CMIP6 prep

---

library(raster)

library(plyr)
library(dplyr)

library(ggplot2)
library(viridisLite)
library(colorspace)
library(RColorBrewer)

library(sf)
library(stringr)

library(ggpubr)
library(kableExtra)

memory.limit(size = 160000)
## [1] Inf

#display.brewer.pal(n = 11, name ="RdYlGn")
colorvec <- brewer.pal(n = 11, name ="RdYlGn")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

1 Data

land_mask <- raster("Masks/land_sea_mask_1degree.nc4") 
land_mask.df <- as.data.frame(land_mask, xy = T) %>% setNames(c("lon","lat","lm")) # 0 if ocean, 1 if land


elev <- raster("Worldclim/wc2.1_10m/wc2.1_10m_elev.tif") 
elev.df <- elev %>% projectRaster(to = land_mask) %>% as.data.frame(xy = T) %>% setNames(c("lon","lat", "z"))


ipcc_regions <- shapefile("Masks/IPCC-WGI-reference-regions-v4.shp") %>% spTransform(crs("EPSG:4326")) 
ipcc_regions.raster <- ipcc_regions %>% rasterize(land_mask)
ipcc_regions.df <- as.data.frame(ipcc_regions.raster, xy = T) %>% setNames(c("lon","lat","Continent","Type","Name","Acronym"))

tas.all_annual <- read.table("CMIP6/tas/tas.all_annual.txt")
pr.all_annual <- read.table("CMIP6/pr/pr.all_annual.txt")
sfcWind.all_annual <- read.table("CMIP6/sfcWind/sfcWind.all_annual.txt")
hfls.all_annual <- read.table("CMIP6/hfls/hfls.all_annual.txt")
hfss.all_annual <- read.table("CMIP6/hfss/hfss.all_annual.txt")
hurs.all_annual <- read.table("CMIP6/hurs/hurs.all_annual.txt")

cmip6_annual <- tas.all_annual %>% merge(pr.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_annual, by = c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_annual, "cmip6_annual.txt")

tas.all_january <- read.table("CMIP6/tas/tas.all_january.txt")
pr.all_january <- read.table("CMIP6/pr/pr.all_january.txt")
sfcWind.all_january <- read.table("CMIP6/sfcWind/sfcWind.all_january.txt")
hfls.all_january <- read.table("CMIP6/hfls/hfls.all_january.txt")
hfss.all_january <- read.table("CMIP6/hfss/hfss.all_january.txt")
hurs.all_january <- read.table("CMIP6/hurs/hurs.all_january.txt")

cmip6_january <- tas.all_january %>% merge(pr.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_january, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_january, "cmip6_january.txt")

tas.all_february <- read.table("CMIP6/tas/tas.all_february.txt")
pr.all_february <- read.table("CMIP6/pr/pr.all_february.txt")
sfcWind.all_february <- read.table("CMIP6/sfcWind/sfcWind.all_february.txt")
hfls.all_february <- read.table("CMIP6/hfls/hfls.all_february.txt")
hfss.all_february <- read.table("CMIP6/hfss/hfss.all_february.txt")
hurs.all_february <- read.table("CMIP6/hurs/hurs.all_february.txt")

cmip6_february <- tas.all_february %>% merge(pr.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_february, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_february, "cmip6_february.txt")

tas.all_march <- read.table("CMIP6/tas/tas.all_march.txt")
pr.all_march <- read.table("CMIP6/pr/pr.all_march.txt")
sfcWind.all_march <- read.table("CMIP6/sfcWind/sfcWind.all_march.txt")
hfls.all_march <- read.table("CMIP6/hfls/hfls.all_march.txt")
hfss.all_march <- read.table("CMIP6/hfss/hfss.all_march.txt")
hurs.all_march <- read.table("CMIP6/hurs/hurs.all_march.txt")

cmip6_march <- tas.all_march %>% merge(pr.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_march, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_march, "cmip6_march.txt")

tas.all_april <- read.table("CMIP6/tas/tas.all_april.txt")
pr.all_april <- read.table("CMIP6/pr/pr.all_april.txt")
sfcWind.all_april <- read.table("CMIP6/sfcWind/sfcWind.all_april.txt")
hfls.all_april <- read.table("CMIP6/hfls/hfls.all_april.txt")
hfss.all_april <- read.table("CMIP6/hfss/hfss.all_april.txt")
hurs.all_april <- read.table("CMIP6/hurs/hurs.all_april.txt")

cmip6_april <- tas.all_april %>% merge(pr.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_april, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_april, "cmip6_april.txt")

tas.all_may <- read.table("CMIP6/tas/tas.all_may.txt")
pr.all_may <- read.table("CMIP6/pr/pr.all_may.txt")
sfcWind.all_may <- read.table("CMIP6/sfcWind/sfcWind.all_may.txt")
hfls.all_may <- read.table("CMIP6/hfls/hfls.all_may.txt")
hfss.all_may <- read.table("CMIP6/hfss/hfss.all_may.txt")
hurs.all_may <- read.table("CMIP6/hurs/hurs.all_may.txt")

cmip6_may <- tas.all_may %>% merge(pr.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_may, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_may, "cmip6_may.txt")

tas.all_june <- read.table("CMIP6/tas/tas.all_june.txt")
pr.all_june <- read.table("CMIP6/pr/pr.all_june.txt")
sfcWind.all_june <- read.table("CMIP6/sfcWind/sfcWind.all_june.txt")
hfls.all_june <- read.table("CMIP6/hfls/hfls.all_june.txt")
hfss.all_june <- read.table("CMIP6/hfss/hfss.all_june.txt")
hurs.all_june <- read.table("CMIP6/hurs/hurs.all_june.txt")

cmip6_june <- tas.all_june %>% merge(pr.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_june, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_june, "cmip6_june.txt")

tas.all_july <- read.table("CMIP6/tas/tas.all_july.txt")
pr.all_july <- read.table("CMIP6/pr/pr.all_july.txt")
sfcWind.all_july <- read.table("CMIP6/sfcWind/sfcWind.all_july.txt")
hfls.all_july <- read.table("CMIP6/hfls/hfls.all_july.txt")
hfss.all_july <- read.table("CMIP6/hfss/hfss.all_july.txt")
hurs.all_july <- read.table("CMIP6/hurs/hurs.all_july.txt")

cmip6_july <- tas.all_july %>% merge(pr.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_july, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_july, "cmip6_july.txt")

tas.all_august <- read.table("CMIP6/tas/tas.all_august.txt")
pr.all_august <- read.table("CMIP6/pr/pr.all_august.txt")
sfcWind.all_august <- read.table("CMIP6/sfcWind/sfcWind.all_august.txt")
hfls.all_august <- read.table("CMIP6/hfls/hfls.all_august.txt")
hfss.all_august <- read.table("CMIP6/hfss/hfss.all_august.txt")
hurs.all_august <- read.table("CMIP6/hurs/hurs.all_august.txt")

cmip6_august <- tas.all_august %>% merge(pr.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_august, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_august, "cmip6_august.txt")

tas.all_september <- read.table("CMIP6/tas/tas.all_september.txt")
pr.all_september <- read.table("CMIP6/pr/pr.all_september.txt")
sfcWind.all_september <- read.table("CMIP6/sfcWind/sfcWind.all_september.txt")
hfls.all_september <- read.table("CMIP6/hfls/hfls.all_september.txt")
hfss.all_september <- read.table("CMIP6/hfss/hfss.all_september.txt")
hurs.all_september <- read.table("CMIP6/hurs/hurs.all_september.txt")

cmip6_september <- tas.all_september %>% merge(pr.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_september, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_september, "cmip6_september.txt")

tas.all_october <- read.table("CMIP6/tas/tas.all_october.txt")
pr.all_october <- read.table("CMIP6/pr/pr.all_october.txt")
sfcWind.all_october <- read.table("CMIP6/sfcWind/sfcWind.all_october.txt")
hfls.all_october <- read.table("CMIP6/hfls/hfls.all_october.txt")
hfss.all_october <- read.table("CMIP6/hfss/hfss.all_october.txt")
hurs.all_october <- read.table("CMIP6/hurs/hurs.all_october.txt")

cmip6_october <- tas.all_october %>% merge(pr.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_october, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_october, "cmip6_october.txt")

tas.all_november <- read.table("CMIP6/tas/tas.all_november.txt")
pr.all_november <- read.table("CMIP6/pr/pr.all_november.txt")
sfcWind.all_november <- read.table("CMIP6/sfcWind/sfcWind.all_november.txt")
hfls.all_november <- read.table("CMIP6/hfls/hfls.all_november.txt")
hfss.all_november <- read.table("CMIP6/hfss/hfss.all_november.txt")
hurs.all_november <- read.table("CMIP6/hurs/hurs.all_november.txt")

cmip6_november <- tas.all_november %>% merge(pr.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_november, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_november, "cmip6_november.txt")

tas.all_december <- read.table("CMIP6/tas/tas.all_december.txt")
pr.all_december <- read.table("CMIP6/pr/pr.all_december.txt")
sfcWind.all_december <- read.table("CMIP6/sfcWind/sfcWind.all_december.txt")
hfls.all_december <- read.table("CMIP6/hfls/hfls.all_december.txt")
hfss.all_december <- read.table("CMIP6/hfss/hfss.all_december.txt")
hurs.all_december <- read.table("CMIP6/hurs/hurs.all_december.txt")

cmip6_december <- tas.all_december %>% merge(pr.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(sfcWind.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfls.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hfss.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>%
  merge(hurs.all_december, by = c("lon","lat","model","period","month","lm", "Continent","Type","Name","Acronym","source")) %>% 
  merge(elev.df, by = c("lon", "lat"))

write.table(cmip6_december, "cmip6_december.txt")

2 Aridity Index

2.1 De Martonne

cmip6 <- read.table("cmip6_annual.txt")

breaks.martonne <- c(0,5,10,20,"30","40",Inf)
cat.martonne <- c("[0,5]" = "Desert", "(5,10]"= "Arid", "(10,20]" = "Semi-arid","(20,30]" = "Temperate", "(30,40]" = "Humid", "(40,Inf]" = "Forest")
col.martonne <- c("Desert" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Temperate" = colorvec[8], "Humid"= colorvec[10],"Forest"= colorvec[11])

cmip6$AIm <- with(cmip6,pr*60*60*24*365/(tas-273.15+10)) # pr in mm/y, tas in Ceslsius


cmip6$cat.AIm <- cmip6$AIm %>% cut(breaks.martonne, include.lowest = T) %>% revalue(cat.martonne)


map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1970_2000" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AIm))+
    scale_fill_manual(values = col.martonne, na.translate = F)+
    labs(title = paste(i, "1970-2000", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

2.2 UNEP

ea from Rh:

ea = 6.108(RH/100)exp(17.27*T/(T+237.3)) Refs https://www.nature.com/articles/s41598-023-40499-6 : Allen and Tetens

id_vars <- c("lon","lat","model","period","lm", "Continent","Type","Name","Acronym","source")

cmip6_annual <- mutate(read.table("cmip6_annual.txt"),
                       t = tas - 273.15,
                       Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                       P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                       es = 0.61078*exp(17.27*t/(t+237.3)),
                       ea = (hurs/100)*es,
                       delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                       gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                       ET0_annual = ((0.408*delta*Rn+gamma*900/((t)+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind)))) 

cmip6_01 <- mutate(read.table("cmip6_january.txt"),
                   t = tas - 273.15,
                   pr_01 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_01 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind)))) # Evapotranspiration, G considered to be negligible. in mm/day. Wind speed at 10 m is converted to wind speed at 2 m by multiplying by 0.748

cmip6_02 <- mutate(read.table("cmip6_february.txt"),
                   t = tas - 273.15,
                   pr_02 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_02 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_03 <- mutate(read.table("cmip6_march.txt"),
                   t = tas - 273.15,
                   pr_03 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_03 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_04 <- mutate(read.table("cmip6_april.txt"),
                   t = tas - 273.15,
                   pr_04 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_04 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_05 <- mutate(read.table("cmip6_may.txt"),
                   t = tas - 273.15,
                   pr_05 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_05 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_06 <- mutate(read.table("cmip6_june.txt"),
                   t = tas - 273.15,
                   pr_06 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_06 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_07 <- mutate(read.table("cmip6_july.txt"),
                   t = tas - 273.15,
                   pr_07 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_07 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_08 <- mutate(read.table("cmip6_august.txt"),
                   t = tas - 273.15,
                   pr_08 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_08 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_09 <- mutate(read.table("cmip6_september.txt"),
                   t = tas - 273.15,
                   pr_09 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_09 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_10 <- mutate(read.table("cmip6_october.txt"),
                   t = tas - 273.15,
                   pr_10 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_10 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_11 <- mutate(read.table("cmip6_november.txt"),
                   t = tas - 273.15,
                   pr_11 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_11 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

cmip6_12 <- mutate(read.table("cmip6_december.txt"),
                   t = tas - 273.15,
                   pr_12 = pr*60*60*24, # mm/day
                   Rn = (hfls + hfss) * 86400 * 1e-6, # conversion from W/m2 to MJ/m2/day.
                   P = 101.3 * ((293-0.0065*z)/293)^(5.26),
                   es = 0.61078*exp(17.27*t/(t+237.3)),
                   ea = (hurs/100)*es,
                   delta = 4098*es/(t+237.3)^2, #slope vapor pressure
                   gamma = 0.00163*P/(2.501-52.6361e-3*t), #psychrometric constant
                   ET0_12 = ((0.408*delta*Rn+gamma*900/(t+273)*0.748*sfcWind*(es-ea))/(delta+gamma*(1+0.34*0.748*sfcWind))))

et0_monthly <- select(cmip6_annual, id_vars) %>%
  merge(select(cmip6_01, c(id_vars, "ET0_01", "pr_01")), by = id_vars) %>%
  merge(select(cmip6_02, c(id_vars, "ET0_02", "pr_02")), by = id_vars) %>%
  merge(select(cmip6_03, c(id_vars, "ET0_03", "pr_03")), by = id_vars) %>%
  merge(select(cmip6_04, c(id_vars, "ET0_04", "pr_04")), by = id_vars) %>%
  merge(select(cmip6_05, c(id_vars, "ET0_05", "pr_05")), by = id_vars) %>%
  merge(select(cmip6_06, c(id_vars, "ET0_06", "pr_06")), by = id_vars) %>%
  merge(select(cmip6_07, c(id_vars, "ET0_07", "pr_07")), by = id_vars) %>%
  merge(select(cmip6_08, c(id_vars, "ET0_08", "pr_08")), by = id_vars) %>%
  merge(select(cmip6_09, c(id_vars, "ET0_09", "pr_09")), by = id_vars) %>%
  merge(select(cmip6_10, c(id_vars, "ET0_10", "pr_10")), by = id_vars) %>%   
  merge(select(cmip6_11, c(id_vars, "ET0_11", "pr_11")), by = id_vars) %>%
  merge(select(cmip6_12, c(id_vars, "ET0_12", "pr_12")), by = id_vars)

et0_monthly$spr <- with(et0_monthly, pr_01*31 + pr_02 * 28 + pr_03 * 31 + pr_04 * 30 + pr_05 * 31 + pr_06 * 30 +
                          pr_07 * 31 + pr_08 * 31 + pr_09 * 30 + pr_10 * 31 + pr_11 * 30 + pr_12 * 31)
et0_monthly$sET0 <- with(et0_monthly, ET0_01*31 + ET0_02 * 28 + ET0_03 * 31 + ET0_04 * 30 + ET0_05 * 31 + ET0_06 * 30 +
                           ET0_07 * 31 + ET0_08 * 31 + ET0_09 * 30 + ET0_10 * 31 + ET0_11 * 30 + ET0_12 * 31)

cmip6 <- merge(cmip6_annual, et0_monthly, by = id_vars)
cmip6$AI <- with(cmip6, abs(spr/sET0))  
cmip6$AI_annual <- with(cmip6, abs((pr*60*60*24*365)/(ET0_annual*365)))

breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")
col.cat <- c("Hyperarid" = colorvec[2], "Arid"= colorvec[4],"Semi-arid"= colorvec[5], "Dry subhumid" = colorvec[8], "Humid"= colorvec[10], "Cold" = "powderblue")

cmip6$cat.AI <- cmip6$AI %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6$cat.AI[which(cmip6$sET0 < 400)] <- "Cold"
cmip6$cat.AI %>% unique()

cmip6$cat.AI <- cmip6$AI %>% cut(breaks = breaks.unesco) %>% 
  revalue(cat.unesco)


cmip6$cat.AI_annual <- cmip6$AI_annual %>% cut(breaks = breaks.unesco) %>% 
  revalue(cat.unesco)

cmip6$cat.AI_annual[which(cmip6$ET0_annual*365 < 400)] <- "Cold"

write.table(cmip6, "cmip6_AI.txt")

3 Maps AI by model, all % calculated without Antarctica

cmip6 <- read.table("cmip6_AI.txt")

breaks.unesco <- c(-Inf,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(-Inf,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid")

3.1 1850-1880

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1850_1880" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "1850-1880", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "1850_1880"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))


tab <- cmip6 %>% subset(period == "1850_1880"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)


knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	5.1	9.0	13.2	6.1	33.4	38.0	28.3	0.3
CESM	5.2	10.8	14.7	5.0	35.7	35.6	28.5	0.3
CMCC	5.6	7.2	9.5	5.3	27.6	42.1	30.1	0.3
CMCC-ESM2	6.0	8.2	7.5	2.8	24.5	24.0	51.2	0.3
CNRM	5.7	12.1	11.4	6.6	35.8	38.0	25.9	0.2
EC-Earth3	8.6	9.0	11.0	4.5	33.1	36.7	30.1	0.2
FGOALS	5.6	12.8	15.9	7.2	41.5	29.9	28.4	0.2
GFDL-ESM4	5.6	9.3	8.7	4.1	27.7	37.0	34.9	0.2
INM	4.4	8.2	13.4	5.8	31.8	40.6	27.4	0.3
INM-CM5	3.5	8.8	11.5	5.6	29.4	43.2	27.1	0.3
MPI	8.5	10.9	10.1	4.4	33.9	34.5	31.4	0.2
MRI	6.3	10.6	8.7	3.8	29.4	39.6	30.9	0.3
NorESM-2-MM	5.1	10.4	18.4	5.9	39.8	35.4	24.4	0.3
mean	5.8	9.8	11.8	5.2	32.6	36.5	30.7	0.3
sd	1.4	1.6	3.2	1.2	7.4	5.1	6.7	0.1

3.2 1970-2000

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1970_2000" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "1850-1880", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

library(ggrepel) 

pie <- cmip6 %>% subset(period == "1970_2000" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  select(c("source", "cat.AI")) %>%
  reshape2::melt(id.vars = c("source")) %>%
  group_by(source, value) %>%
  summarise(count = n()) %>%
  mutate(perc = round(count/sum(count)*100, 1),
         lab.y = rev(cumsum(rev(count))) - count*0.5) %>% 
  ungroup() 

pie_list <- list()

for(i in unique(pie$source)){
  p <- ggplot(subset(pie, source == i))+
    geom_bar(aes(x = "", y = count, fill = value), width = 1, stat = "identity")+
    coord_polar("y", start = 0)+
    scale_fill_manual(values = col.cat)+
    geom_label_repel(aes(x = "", y = lab.y, label = perc), nudge_x = 0.5)+
    labs(fill = "", title = i)+
    theme_void()+theme(legend.position = "right")
  pie_list[[i]] <- p
}

ggpubr::ggarrange(plotlist = pie_list, nrow = 4, ncol = 4, common.legend = T, legend = "bottom")

pie <- cmip6 %>% subset(period == "1970_2000" & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  select(c("cat.AI")) %>%
  table() %>% as.data.frame() %>%
  mutate(perc = round(Freq/sum(Freq)*100, 1),
         lab.y = rev(cumsum(rev(Freq))) - Freq*0.5) 

p <- ggplot(pie)+
    geom_bar(aes(x = "", y = Freq, fill = cat.AI), width = 1, stat = "identity")+
    coord_polar("y", start = 0)+
    scale_fill_manual(values = col.cat)+
    geom_label_repel(aes(x = "", y = lab.y, label = perc), nudge_x = 0.5)+
    labs(fill = "", title = "Multimodel CMIP6")+
    theme_void()+theme(legend.position = "right")

print(p)

sdmm <- cmip6 %>% subset(period == "1970_2000"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))


tab <- cmip6 %>% subset(period == "1970_2000"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)


knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	5.2	8.6	13.2	6.3	33.3	37.4	29.1	0.3
CESM	5.3	10.7	14.1	5.2	35.3	35.0	29.5	0.3
CMCC	5.1	7.4	10.1	5.4	28.0	42.2	29.5	0.3
CMCC-ESM2	6.4	7.9	7.7	2.5	24.5	24.4	50.8	0.3
CNRM	5.8	12.2	11.3	6.5	35.8	38.3	25.7	0.2
EC-Earth3	8.4	9.5	10.8	4.4	33.1	35.7	31.0	0.2
FGOALS	5.9	10.9	10.0	4.9	31.7	25.1	42.9	0.2
GFDL-ESM4	5.8	8.9	8.9	4.1	27.7	35.3	36.8	0.2
INM	4.4	8.4	13.3	5.7	31.8	40.1	27.8	0.3
INM-CM5	3.4	8.5	12.0	5.6	29.5	43.0	27.3	0.3
MPI	8.7	11.1	10.2	4.3	34.3	34.0	31.5	0.2
MRI	6.1	10.8	9.2	3.9	30.0	38.9	30.9	0.3
NorESM-2-MM	5.2	10.1	17.1	5.6	38.0	33.4	28.3	0.3
mean	5.8	9.6	11.4	5.0	31.8	35.6	32.4	0.3
sd	1.4	1.5	2.5	1.1	6.5	5.7	7.1	0.1

3.3 1985-2010

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "1985_2015" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "1985-2010", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "1985_2015"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))


tab <- cmip6 %>% subset(period == "1985_2015"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)


knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	5.3	8.4	13.8	6.6	34.1	37.2	28.3	0.3
CESM	4.8	10.8	14.8	5.4	35.8	35.4	28.5	0.3
CMCC	5.2	7.6	10.4	5.8	29.0	42.9	27.9	0.3
CMCC-ESM2	6.2	8.5	7.7	2.5	24.9	24.5	50.3	0.3
CNRM	5.8	12.1	11.3	6.2	35.4	39.2	25.2	0.2
EC-Earth3	8.3	9.8	10.7	4.4	33.2	36.7	29.9	0.2
FGOALS	5.8	11.4	10.4	4.9	32.5	25.3	42.0	0.2
GFDL-ESM4	5.9	9.1	9.1	3.9	28.0	36.5	35.3	0.2
INM	4.0	7.7	11.3	5.2	28.2	36.5	35.0	0.3
INM-CM5	3.5	8.3	12.0	5.4	29.2	40.6	29.9	0.3
MPI	8.7	11.1	10.3	4.4	34.5	34.1	31.0	0.2
MRI	6.2	10.7	9.1	3.8	29.8	38.4	31.6	0.3
NorESM-2-MM	4.7	9.9	15.9	5.6	36.1	30.3	33.2	0.3
mean	5.7	9.6	11.3	4.9	31.5	35.2	32.9	0.3
sd	1.5	1.5	2.3	1.1	6.4	5.5	6.7	0.1

3.4 SSP 2-4.5

3.4.1 2030-2060

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP245" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2030_2060" & model == "SSP245" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2030_2060" & model == "SSP245"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	5.0	9.2	14.6	6.5	35.3	38.1	26.2	0.3
CESM	4.7	12.4	15.3	5.6	38.0	35.6	26.2	0.3
CMCC	4.4	8.4	11.0	6.1	29.9	44.3	25.5	0.3
CMCC-ESM2	6.1	9.1	8.0	2.8	26.0	25.3	48.4	0.3
CNRM	6.0	12.7	12.1	6.8	37.6	39.8	22.3	0.2
EC-Earth3	7.5	9.2	11.1	4.5	32.3	36.3	31.2	0.2
FGOALS	6.9	13.3	18.5	7.8	46.5	30.3	22.9	0.2
GFDL-ESM4	5.4	10.1	9.5	4.4	29.4	36.9	33.4	0.2
INM	4.6	8.4	13.8	5.9	32.7	40.4	26.7	0.3
INM-CM5	3.5	8.6	12.2	5.4	29.7	43.1	26.9	0.3
MPI	9.0	11.9	11.1	4.5	36.5	35.1	28.2	0.2
MRI	6.2	11.2	9.3	4.1	30.8	40.0	28.9	0.3
NorESM-2-MM	5.3	11.8	16.9	5.6	39.6	33.5	26.7	0.3
mean	4.0	7.3	8.7	3.7	23.7	25.6	48.0	2.6
sd	1.0	1.2	2.2	0.9	5.3	3.6	4.6	0.0

3.4.2 2070-2100

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP245" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2070_2100" & model == "SSP245" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2070_2100" & model == "SSP245"   & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	4.8	10.5	15.0	7.3	37.6	38.2	24.1	0.3
CESM	6.1	11.6	15.4	5.4	38.5	36.9	24.3	0.3
CMCC	4.3	8.0	12.3	6.4	31.0	46.6	22.1	0.3
CMCC-ESM2	6.2	9.5	8.5	2.8	27.0	26.8	46.0	0.3
CNRM	6.3	12.9	12.4	6.8	38.4	41.0	20.4	0.2
EC-Earth3	7.6	9.1	11.0	4.4	32.1	37.3	30.3	0.2
FGOALS	6.9	14.0	18.1	8.1	47.1	30.3	22.3	0.2
GFDL-ESM4	5.9	9.7	9.8	4.4	29.8	38.4	31.6	0.2
INM	4.8	8.2	14.1	5.8	32.9	40.6	26.2	0.3
INM-CM5	3.7	9.5	11.2	5.9	30.3	43.5	25.9	0.3
MPI	9.0	11.9	11.1	4.5	36.5	35.1	28.2	0.2
MRI	6.5	11.0	9.6	4.0	31.1	40.8	27.9	0.3
NorESM-2-MM	5.9	11.7	17.0	5.9	40.5	34.0	25.3	0.3
mean	4.2	7.4	8.8	3.8	24.2	26.2	47.0	2.6
sd	1.0	1.3	2.1	1.0	5.4	3.7	4.5	0.0

3.5 SSP 3-7.0

3.5.1 2030-2060

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP370" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2030_2060" & model == "SSP370" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2030_2060" & model == "SSP370"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	4.9	9.6	14.3	6.8	35.6	37.4	26.8	0.3
CESM	5.2	11.4	15.6	5.4	37.6	35.7	26.5	0.3
CMCC	4.2	8.3	11.3	5.9	29.7	44.1	25.9	0.3
CMCC-ESM2	3.3	8.4	9.0	3.6	24.3	26.9	48.6	0.3
CNRM	6.2	12.9	12.9	7.1	39.1	39.6	21.0	0.2
EC-Earth3	7.9	9.7	11.6	4.8	34.0	38.5	27.3	0.2
FGOALS	6.1	11.9	11.9	5.3	35.2	27.5	37.0	0.2
GFDL-ESM4	6.3	9.4	9.9	4.1	29.7	35.9	34.1	0.2
INM	4.6	8.3	13.5	5.9	32.3	40.3	27.1	0.3
INM-CM5	3.4	8.8	12.1	5.6	29.9	43.2	26.5	0.3
MPI	8.7	12.2	11.0	4.5	36.4	36.4	26.9	0.2
MRI	6.5	10.8	9.4	3.6	30.3	40.2	29.2	0.3
NorESM-2-MM	5.4	11.0	17.9	5.7	40.0	33.2	26.5	0.3
mean	3.9	7.1	8.6	3.7	23.3	25.6	48.5	2.6
sd	1.1	1.1	1.8	0.8	4.8	3.6	4.8	0.0

3.5.2 2070-2100

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP370" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2070_2100" & model == "SSP370" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2070_2100" & model == "SSP370"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	5.5	10.5	15.4	6.9	38.3	37.3	24.2	0.3
CESM	6.0	12.2	15.6	5.7	39.5	37.4	22.9	0.3
CMCC	5.2	9.1	12.5	7.3	34.1	44.1	21.6	0.3
CMCC-ESM2	6.9	9.9	8.1	3.0	27.9	25.9	45.9	0.3
CNRM	7.3	12.9	13.2	6.7	40.1	39.3	20.3	0.2
EC-Earth3	7.6	9.8	12.4	4.7	34.5	40.3	24.9	0.2
FGOALS	6.7	13.3	13.3	4.4	37.7	28.0	34.1	0.2
GFDL-ESM4	7.5	8.5	10.9	4.7	31.6	37.0	31.1	0.2
INM	6.0	8.2	14.7	6.1	35.0	39.2	25.5	0.3
INM-CM5	3.9	10.3	14.1	6.1	34.4	40.4	24.9	0.3
MPI	8.7	12.2	11.0	4.5	36.4	36.4	26.9	0.2
MRI	7.2	11.9	9.9	4.1	33.1	39.0	27.6	0.3
NorESM-2-MM	5.9	11.8	17.5	5.9	41.1	33.7	24.9	0.3
mean	4.5	7.5	9.0	3.8	24.8	25.6	47.0	2.6
sd	0.9	1.2	1.8	0.9	4.8	3.5	4.7	0.0

3.6 SSP 5-8.5

3.6.1 2030-2060

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2030_2060" & model == "SSP585" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2030-2060", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2030_2060" & model == "SSP585" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2030_2060" & model == "SSP585"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(14,15), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	4.8	9.2	15.0	7.5	36.5	37.6	25.6	0.3
CESM	5.8	11.5	15.6	5.2	38.1	36.0	25.7	0.3
CMCC	4.4	8.0	11.5	6.5	30.4	44.8	24.5	0.3
CMCC-ESM2	6.5	9.3	8.0	2.5	26.3	25.4	48.0	0.3
CNRM	6.2	13.3	12.4	6.9	38.8	40.0	21.0	0.2
EC-Earth3	7.8	9.4	11.8	4.6	33.6	38.8	27.4	0.2
FGOALS	6.8	13.5	18.1	7.7	46.1	30.8	22.8	0.2
GFDL-ESM4	6.1	9.4	9.7	4.7	29.9	36.4	33.4	0.2
INM	4.6	8.4	14.0	6.0	33.0	40.6	26.2	0.3
INM-CM5	3.5	8.6	12.4	5.8	30.3	43.6	25.8	0.3
MPI	8.6	12.4	11.1	4.5	36.6	37.0	26.2	0.2
MRI	6.0	11.1	9.6	4.0	30.7	40.8	28.1	0.3
NorESM-2-MM	5.5	11.4	17.6	5.9	40.4	33.2	26.2	0.3
mean	4.1	7.3	8.9	3.8	24.1	26.0	47.3	2.6
sd	1.0	1.3	2.1	1.0	5.4	3.7	4.7	0.0

3.6.2 2070-2100

map_list <- list()

for(i in unique(cmip6$source)){
  g <- ggplot(data = subset(cmip6, period == "2070_2100" & model == "SSP585" & source == i)) + geom_raster(aes(x=lon, y = lat,  fill = cat.AI))+
    scale_fill_manual(values = col.cat, na.translate = F)+
    labs(title = paste(i, "2070-2100", sep = ", "), fill = "")+
    theme_void()
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 3, nrow = 5, common.legend = T, legend = "bottom")

sdmm <- cmip6 %>% subset(period == "2070_2100" & model == "SSP585" &!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  ungroup() %>% group_by(cat.AI) %>% summarise(mean = round(mean(percent, na.rm = T),1), sd = round(sd(percent, na.rm = T),1)) %>%
  t() %>% as.data.frame() %>% setNames(c(.[1,1:6], "NA")) %>% .[-1, ] %>% sapply(as.numeric) %>% as.data.frame() %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")]), "source" = c("mean","sd")) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

tab <- cmip6 %>% subset(period == "2070_2100" & model == "SSP585"  & !Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55) %>%
  group_by(source, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(source) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(source~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("source","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA")) %>% 
  rbind(sdmm)

knitr::kable(tab, digits = 2, escape = F)  %>%
  kable_styling(bootstrap_options = "bordered") %>% 
  column_spec(c(6,9), italic = T, include_thead = T) %>% row_spec(c(13,14), bold = T) 

source	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
CAS-ESM2	6.1	10.8	15.9	6.5	39.3	37.5	22.9	0.3
CMCC	4.4	9.3	13.5	6.7	33.9	46.7	19.2	0.3
CMCC-ESM2	5.6	10.5	8.9	2.9	27.9	28.1	43.8	0.3
CNRM	6.0	13.6	12.8	6.4	38.8	40.9	20.1	0.2
EC-Earth3	7.7	9.6	13.0	4.9	35.2	40.5	23.9	0.2
FGOALS	15.6	15.0	18.3	6.8	55.7	24.4	19.8	0.2
GFDL-ESM4	6.5	10.0	10.8	4.8	32.1	37.6	30.1	0.2
INM	4.4	8.7	14.8	6.0	33.9	41.7	24.1	0.3
INM-CM5	3.9	9.8	12.0	6.0	31.7	44.1	23.9	0.3
MPI	8.7	12.3	11.1	4.7	36.8	36.8	26.2	0.2
MRI	7.1	11.5	10.6	4.1	33.3	40.1	26.4	0.3
NorESM-2-MM	6.4	12.8	17.4	6.5	43.1	34.3	22.3	0.3
mean	4.8	7.8	9.2	3.8	25.6	26.2	45.5	2.6
sd	2.1	1.3	2.0	0.8	6.2	4.4	4.6	0.0

4 Agreement between models

4.1 Summary dataframe

4.1.1 Difference with AI ref for all SSP

test <- cmip6ext %>% subset(model %in% c("historical","SSP585") & period %in% c("1970_2000", "2070_2100") & lon == 25.5 & lat == 10.5) %>% select(c("source","lon","lat","model","period","AI","cat.AI","AI.ref"))
                                                                               
                                                                               
test2 <- cmip6s %>% subset(model == "SSP585" & period == "2070_2100" & lon == 25.5 & lat == 10.5) %>% 
  select(c("model","period","AI.mean","cat.AI","AI.ref.mean","diff.AI.mean", "nb.models", "diff.AI.mean"))

test2

cmip6 <- read.table("cmip6_AI.txt")

cmip6.ref <- cmip6 %>% filter(period == "1970_2000") %>% select(c("lon","lat","model","period","source","AI","spr","t")) %>% setNames(c("lon","lat","model","period","source","AI.ref","spr.ref","t.ref"))

cmip6ext <- merge(cmip6, select(cmip6.ref, c("lon","lat","source","AI.ref","spr.ref","t.ref")) , by = c("lon","lat","source"), all = T) %>%  mutate(diff.AI = AI - AI.ref, diff.pr = spr - spr.ref, diff.t = t - t.ref)

write.table(cmip6ext, "cmip6ext.txt")

cmip6 <- read.table("cmip6ext.txt")
cmip6s <- cmip6 %>% 
  group_by(lon,lat,Continent, Type, Name, Acronym, lm, model, period, z) %>% 
  dplyr::summarise(t.mean = mean(t, na.rm = T), t.sd = sd(t, na.rm = T), # t in °C
                   pr.mean = mean(spr, na.rm = T) , pr.sd = sd(spr, na.rm = T),  # pr in mm/y
                   sfcWind.mean = mean(sfcWind, na.rm = T), sfcWind.sd = sd(sfcWind, na.rm = T),
                   Rn.mean = mean(Rn, na.rm = T), Rn.sd = sd(Rn, na.rm = T), # in MJ/m2/day
                   ea.mean = mean(ea, na.rm = T), ea.sd = sd(ea, na.rm = T), # in kPa
                   ET0 = mean(sET0, na.rm = T), ET0.sd = sd(sET0, na.rm = T), # in mm/y
                   AI.mean = mean(AI, na.rm = T), AI.sd = sd(AI, na.rm = T),
                   AI.ref.mean = mean(AI.ref, na.rm = T), AI.ref.sd = sd(AI.ref, na.rm = T),
                   diff.AI.mean = mean(diff.AI, na.rm = T), diff.AI.sd = sd(diff.AI, na.rm = T),
                   spr.ref.mean = mean(spr.ref, na.rm = T), spr.ref.sd = sd(spr.ref, na.rm = T),
                   diff.pr.mean = mean(diff.pr, na.rm = T), diff.pr.sd = sd(diff.pr, na.rm = T),
                   t.ref.mean = mean(t.ref, na.rm = T), t.ref.sd = sd(t.ref, na.rm =T),
                   diff.t.mean = mean(diff.t, na.rm = T), diff.t.sd = sd(diff.t, na.rm = T),
                   Hyperarid = sum(cat.AI == "Hyperarid"),
                   Arid = sum(cat.AI == "Arid"),
                   'Semi-arid' = sum(cat.AI == "Semi-arid"),
                   'Dry subhumid' = sum(cat.AI == "Dry subhumid") ,
                   'Humid' = sum(cat.AI == "Humid"),
                   'Cold' = sum(cat.AI == "Cold"),
                   'NA' = sum(is.na(cat.AI)),
                   nb.models = n(),
                    plus.AI = table(diff.AI > 0)["TRUE"],
                    minus.AI = table(diff.AI < 0)["TRUE"]) %>%
ungroup() %>% as.data.frame()

cmip6s$cat.maj <- colnames(select(cmip6s, c("Hyperarid","Arid",'Semi-arid',"Dry subhumid",'Humid',"Cold","NA")))[max.col(select(cmip6s, c("Hyperarid","Arid",'Semi-arid',"Dry subhumid",'Humid',"Cold","NA")),ties.method = "first")]
cmip6s$nb.maj <- cmip6s %>% select(c("Hyperarid","Arid",'Semi-arid',"Dry subhumid",'Humid',"Cold","NA")) %>% apply(1, max)
cmip6s$prop.maj <- with(cmip6s, nb.maj / nb.models)



breaks.unesco <- c(-Inf,0,0.03,0.2,0.5,0.65,Inf)
cat.unesco <- c("(0,0.03]"= "Hyperarid", "(0.03,0.2]" = "Arid", "(0.2,0.5]" = "Semi-arid", "(0.5,0.65]" = "Dry subhumid","(0.65, Inf]" = "Humid", "(-Inf,0]" = "Cold")

cmip6s$cat.AI <- cmip6s$AI.mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6s$cat.AI[which(cmip6s$ET0 < 400)] <- "Cold"

cmip6s$cat.AI.ref <- cmip6s$AI.ref.mean %>% cut(breaks = breaks.unesco) %>% revalue(cat.unesco)
cmip6s$cat.AI.ref[which(cmip6s$ET0 < 400)] <- "Cold"

cmip6s$same.catAI <- with(cmip6s, ifelse(cat.maj == cat.AI, "same", "diff"))

cmip6s$cat.dryer <- with(cmip6s, 
          ifelse(cat.AI.ref == "Arid" & cat.AI %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter", 
          ifelse(cat.AI.ref == "Arid" & cat.AI == "Hyperarid", "dryer",
          ifelse(cat.AI.ref == "Semi-arid" & cat.AI %in% c("Cold","Humid","Dry subhumid"), "wetter",
          ifelse(cat.AI.ref == "Semi-arid" & cat.AI %in% c("Hyperarid","Arid"), "dryer",
          ifelse(cat.AI.ref == "Dry subhumid" & cat.AI %in% c("Cold","Humid"), "wetter",
          ifelse(cat.AI.ref == "Dry subhumid" & cat.AI %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.ref == "Humid" & cat.AI == "Cold","wetter",
          ifelse(cat.AI.ref == "Humid" & cat.AI %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.ref == "Cold" & cat.AI %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
          ifelse(cat.AI.ref == "Hyperarid" & cat.AI %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
          ifelse(cat.AI.ref == cat.AI, "same",
          NA))))))))))))

cmip6sref <- cmip6s %>% filter(period == "1970_2000") %>% select(c("lon","lat","model","period","cat.maj")) %>% setNames(c("lon","lat","model","period","cat.maj.ref"))

cmip6s <- merge(cmip6s, select(cmip6sref, c("lon","lat", "model","period","cat.maj.ref")), 
                 by = c("lon","lat","model","period"), all = T) 

cmip6s$cat.maj.dryer <- with(cmip6s, 
        ifelse(cat.maj.ref == "Arid" & cat.maj %in% c("Cold","Humid","Dry subhumid","Semi-arid"), "wetter", 
        ifelse(cat.maj.ref == "Arid" & cat.maj == "Hyperarid", "dryer",
        ifelse(cat.maj.ref == "Semi-arid" & cat.maj %in% c("Cold","Humid","Dry subhumid"), "wetter",
        ifelse(cat.maj.ref == "Semi-arid" & cat.maj %in% c("Hyperarid","Arid"), "dryer",
        ifelse(cat.maj.ref == "Dry subhumid" & cat.maj %in% c("Cold","Humid"), "wetter",
        ifelse(cat.maj.ref == "Dry subhumid" & cat.maj %in% c("Semi-arid","Arid","Hyperarid"), "dryer",
        ifelse(cat.maj.ref == "Humid" & cat.maj == "Cold","wetter",
        ifelse(cat.maj.ref == "Humid" & cat.maj %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
        ifelse(cat.maj.ref == "Cold" & cat.maj %in% c("Humid","Dry subhumid","Semi-arid","Arid","Hyperarid"), "dryer",
        ifelse(cat.maj.ref == "Hyperarid" & cat.maj %in% c("Cold","Humid","Dry subhumid","Semi-arid","Arid"), "wetter",
        ifelse(cat.maj.ref == cat.maj, "same",
        NA))))))))))))


write.table(cmip6s, "cmip6s.txt")

cmip6s <- read.table("cmip6s.txt")

cmip6s$Continent[which(cmip6s$Continent == "Europe-Africa")] <- "Mediterranea"

4.2 Compare cat mean and cat maj

cat mean = category of aridity obtained by using the multimodel mean AI cat maj = category of aridity obtained by using the category majoritarily obtained in the 13 models

4.2.1 Where are cat mean and cat maj different

cmip6s %>% group_by(model, period) %>% summarise(maj.equal.mean = table(same.catAI)[2], maj.diff.mean = table(same.catAI)[1], prop.same = round(maj.equal.mean/n()*100,1)) %>% knitr::kable(col.names = c("Model","Period","Number of gridcells in which cat.maj = cat.mean", "Number of gridcells in which cat.maj /= cat.mean", "Proportion of cells with the same category (%)")) %>% kableExtra::kable_styling(bootstrap_options = "bordered")

Model	Period	Number of gridcells in which cat.maj = cat.mean	Number of gridcells in which cat.maj /= cat.mean	Proportion of cells with the same category (%)
SSP245	2030_2060	19884	1798	89.3
SSP245	2070_2100	19890	1792	89.3
SSP370	2030_2060	19988	1694	89.8
SSP370	2070_2100	19911	1771	89.4
SSP585	2030_2060	19799	1883	88.9
SSP585	2070_2100	19796	1886	88.9
historical	1850_1880	19907	1775	89.4
historical	1970_2000	19880	1802	89.3
historical	1985_2015	19936	1746	89.6

cmip6s$cat.AI <- factor(cmip6s$cat.AI, levels =  c("Hyperarid","Arid","Semi-arid","Dry subhumid", "Humid","Cold"))
cmip6s$cat.maj <- factor(cmip6s$cat.maj, levels =  c("Hyperarid","Arid","Semi-arid","Dry subhumid", "Humid","Cold"))

ggplot(subset(cmip6s, cat.AI != "Cold"))+
  geom_col(aes(x= period, y = cat.AI, fill = cat.AI), just = -0.2, width = 0.3)+
  geom_col(aes(x = period, y = cat.maj, fill = cat.maj), just = 1.2, width = 0.3)+
  scale_fill_manual(values = col.cat, na.translate = F)+
  labs(fill = "", y ="")+
  facet_grid(rows = vars(model), switch = "y")+
  theme_minimal()+theme(axis.text.y = element_blank())

tab <- rbind(table(cmip6s$cat.AI)/length(cmip6s$cat.AI)*100, table(cmip6s$cat.maj)/length(cmip6s$cat.maj)*100) %>% t() %>% as.data.frame() %>% setNames(c("Cat mean", "Cat maj"))
knitr::kable(tab, digits =2, caption = "Proportion of gridcells in each aridity category") %>% kable_styling(bootstrap_options = "bordered")

Proportion of gridcells in each aridity category

	Cat mean	Cat maj
Hyperarid	3.51	4.18
Arid	6.25	7.43
Semi-arid	7.63	9.19
Dry subhumid	3.39	0.84
Humid	28.75	28.53
Cold	47.90	47.23

map_list <- list()

for(i in unique(cmip6s$model)){
  df <- subset(cmip6s, model == i)
  for(j in unique(df$period)){
    index <- paste(i, j, sep = "")
    g <- ggplot()+geom_raster(data = subset(df, period == j & same.catAI == "diff"), aes(x = lon, y = lat, fill = cat.AI))+
      borders(colour = "grey60")+
      scale_fill_manual(values =  col.cat, na.translate = F)+
      labs(title = paste(i, j , sep = " "), fill = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }
}

col1 <- ggarrange(plotlist = map_list, nrow = 9, ncol = 1, common.legend = T, legend = "bottom") %>% annotate_figure(top = "cat.mean for gridcells in which\ncat.maj and cat.mean are different")

map_list2 <- list()

for(i in unique(cmip6s$model)){
  df <- subset(cmip6s, model == i)
  for(j in unique(df$period)){
    index <- paste(i, j, sep = "")
    g <- ggplot()+geom_raster(data = subset(df, period == j & same.catAI == "diff"), aes(x = lon, y = lat, fill = cat.maj))+
      borders(colour = "grey60")+
      scale_fill_manual(values =  col.cat, na.translate = F)+
      labs(title = paste(i, j , sep = " "), fill = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list2[[index]] <- g
  }
}

col2 <- ggarrange(plotlist = map_list2, nrow = 9, ncol = 1, common.legend = T, legend = "bottom") %>% annotate_figure(top = "cat.maj for gridcells in which\ncat.maj and cat.mean are different")

ggarrange(col1,col2,ncol = 2, common.legend = T)

4.2.2 Plot cat maj and cat mean for reference period

ggpubr::ggarrange(plotlist = list(
  
  ggplot() + 
    geom_raster(data = subset(cmip6s, period == "1970_2000"), aes(x=lon, y = lat,  fill = cat.maj))+
    geom_point(data = subset(cmip6s, period == "1970_2000" & nb.maj > 10), aes(x =lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
    scale_fill_manual(values = col.cat)+
    labs(title = "Category chosen in majority by models", fill = "")+
    ylim(-55,90)+
    theme_void(base_size = 15)+
    theme(legend.position = "none"),
  
  ggplot() + 
    geom_raster(data = subset(cmip6s, period == "1970_2000"), aes(x=lon, y = lat,  fill = cat.AI))+
   geom_point(data = subset(cmip6s, period == "1970_2000" & nb.maj > 10), aes(x =lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
    scale_fill_manual(values = col.cat)+
    labs(title = "Multimodel average categories of AI", fill = "")+
    ylim(-55,90)+
    theme_void(base_size = 15)+
    theme(legend.position = "none")),
  
  ncol = 2)

map_list <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = cat.maj))+
    geom_point(data = subset(cmip6s, period == i & model == "historical" & nb.maj > 10), aes(x= lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list[[i]] <- g
}


map_list2 <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI))+
   # geom_point(data = subset(cmip6s, period == i & model == "historical" & nb.maj > 10), aes(x= lon, y = lat), size = 0.5, col = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list2[[i]] <- g
}


ggarrange(
  annotate_figure(ggarrange(plotlist = map_list, nrow = 3, common.legend = T, legend = "bottom"), top = "Majority category among models"),
  annotate_figure(ggarrange(plotlist = map_list2, nrow = 3, common.legend = T, legend = "bottom"), top = "Multimodel average category"),
  ncol = 2, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = cat.maj))+
         geom_point(data = subset(cmip6s, period == j & model == i & nb.maj > 10), aes(x= lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
      scale_fill_manual(values = col.cat, na.translate = F)+
      labs(fill = "Cat AI maj", title = paste(i, j, sep = ", "))+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}


map_list2 <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = cat.AI))+
    #     geom_point(data = subset(cmip6s, period == j & model == i & nb.maj > 10), aes(x= lon, y = lat), shape = "'", size = 0.5, col = "grey60")+
      scale_fill_manual(values = col.cat, na.translate = F)+
      labs(fill = "Cat AI mean", title = paste(i, j, sep = ", "))+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list2[[index]] <- g
  }}

ggarrange(plotlist = map_list, nrow = 3, ncol = 2, common.legend = T, legend = "bottom")

ggarrange(plotlist = map_list2, nrow = 3, ncol = 2, common.legend = T, legend = "bottom")

4.3 Multimodel maps of cat mean

cmip6s$cat.AI <- factor(cmip6s$cat.AI, levels = c("Hyperarid","Arid","Semi-arid","Dry subhumid","Humid","Cold"))

map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + 
    geom_tile(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = cat.AI))+
    geom_point(data = subset(cmip6s,nb.maj > 10), aes(x = lon, y = lat), shape = ".", col = "grey20", size = 0.1)+
    borders(colour = "grey60")+
    scale_fill_manual(values =  col.cat, na.translate = F)+
    labs(title = paste("Category AI",i , sep = " "), fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + 
      geom_tile(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = cat.AI))+
      geom_point(data = subset(cmip6s,nb.maj > 10), aes(x = lon, y = lat), shape = ".", col = "grey20", size = 0.1)+
      borders(colour = "grey60")+
      scale_fill_manual(values =  col.cat, na.translate = F)+
      labs(title = paste("Category AI",index, sep = " "), fill = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

5 Maps multimodel AI

quantile(cmip6s$AI.mean, probs = seq(0,1,0.1), na.rm = T) 
##           0%          10%          20%          30%          40%          50% 
## 2.902197e-03 1.993698e-01 5.922059e-01 1.019621e+00 1.457766e+00 2.127218e+00 
##          60%          70%          80%          90%         100% 
## 3.578980e+00 1.168470e+01 2.438867e+01 4.856852e+01 5.329875e+06

ai.breaks <- c(0,0.03,0.2,0.5,0.65,1,10,20,30)
colscale <- colorvec[-c(5,6)]

map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = AI.mean))+
    borders(colour = "grey60")+
    binned_scale(aesthetics = "fill", breaks = ai.breaks, palette = function(x) colscale,
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "Aridity index")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = AI.mean))+
      borders(colour = "grey60")+
      binned_scale(aesthetics = "fill", breaks = ai.breaks, palette = function(x) colscale,
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Aridity index", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

6 AI changes over time and scenarios

v<-NULL
v2<-NULL
for (N in 5:30) {
  v<-c(v,sum(sapply(floor(2*N/3):N,function(i) choose(N,i)*(1/2)^N*2)))# proba d'avoir 2/3 de même signe à partir de N modèles, lorsque N change
  v2<-c(v2,sum(sapply(floor(2*N/3):N,function(i) choose(N,i)*(1/2)^N)))# proba d'avoir 2/3 de signe >0 à partir de N modèles, lorsque N change
}


plot(5:30,v,type="l",ylab="Proba(>2/3 models...)",xlab="Nb of models")
lines(5:30,v2,col="red")

6.1 Changes AI for summary CMIP6

6.2 Boxplot of changes AI

ggpubr::ggarrange(plotlist = list(
  ggplot(subset(cmip6s, diff.AI.mean < 1 & diff.AI.mean > -1 & model == "historical"))+
    geom_boxplot(aes(x=period, y = - diff.AI.mean, col = period), outliers = F)+
    scale_color_manual(values = c("#1A9850","#66BD63","#A6D96A"))+
    labs(x = "", y = "Difference of AI between each period and\nthe reference period 1970-2000")+
    facet_grid(cols = vars(model))+
    theme_minimal() + theme(legend.position = "none", axis.text.x = element_text(angle = 45, vjust = 1, hjust=1)),
  ggplot(subset(cmip6s, diff.AI.mean < 1 & diff.AI.mean > -1 & model != "historical"))+
    geom_boxplot(aes(x=period, y = diff.AI.mean, col = period), outliers = F)+
    scale_color_manual(values = c("#FDAE61","#F46D43"))+
    scale_y_continuous(position = "right")+
    labs(y = "", x = "")+
    facet_grid(cols = vars(model))+
    theme_minimal() + theme(legend.position = "none", axis.text.x = element_text(angle = 45, vjust = 1, hjust=1))),
  widths = c(1.5,3),
  ncol = 2, nrow = 1
)

6.3 Strictly wetter or dryer (reference period: 1970-2000)

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
  index <- paste(i, j, sep = ", ")
  g <- ggplot()+
  geom_point(data = subset(cmip6s, model == i & period == j & plus.AI > 10), aes(x = lon, y = lat, col = "plus"), shape = "+")+
  geom_point(data = subset(cmip6s, model == i & period == j & minus.AI > 10), aes(x = lon, y = lat, col = "minus"), shape = "-")+
    scale_color_manual(values = colorvec[c(1,11)])+
    labs(color = "", x = "", y = "", title = index)+
  theme_minimal()
  map_list[[index]] <- g
  }
}

ggarrange(plotlist = map_list, ncol = 2, nrow = 3, common.legend = T)

map_list <- list()

for(i in c("1850_1880","1985_2015")){
    index <- paste(i, j, sep = " , ")
  g <- ggplot() + geom_tile(data = subset(cmip6s, period == i & model == "historical" & diff.AI.mean > 0), aes(x=lon, y = lat,  fill = "wetter"))+
    geom_tile(data = subset(cmip6s, period == i & model == "historical" & diff.AI.mean < 0), aes(x=lon, y = lat,  fill = "dryer"))+
    geom_point(data = subset(cmip6s, period == i & model == "historical" & plus.AI > 10), aes(x=lon, y = lat), shape = "+", size = 0.5)+
     geom_point(data = subset(cmip6s, period == i & model == "historical" & minus.AI > 10), aes(x=lon, y = lat), shape = "-", size = 0.5)+
    borders(colour = "grey60")+
    scale_fill_manual(values =  c("#F46D43", "#74ADD1"), na.translate = F)+
    labs(title = paste("AI changes between",i , "and 1970-2000", sep = " "), fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")
  map_list[[i]] <- g
  }

ggarrange(plotlist = map_list, ncol = 2)

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & diff.AI.mean > 0), aes(x=lon, y = lat,  fill = "wetter"))+
      geom_tile(data = subset(cmip6s, period == j & model == i & diff.AI.mean < 0), aes(x=lon, y = lat,  fill = "dryer"))+
       geom_point(data = subset(cmip6s, period == j & model == i & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s, period == j & model == i & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
      borders(colour = "grey60")+
      scale_fill_manual(values = c("#F46D43", "#74ADD1"), na.translate = F)+
      labs(fill = "Compared to 1970-2000", title = i)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}



ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

6.4 Wetter or dryer, binned (reference period: 1970-2000)

ano.quantiles <- cmip6s %>% group_by(model, period) %>% summarise(q10 = quantile(diff.AI.mean, probs = 0.1, na.rm =T), 
                                                                  q25 = quantile(diff.AI.mean, probs = 0.25, na.rm =T), 
                                                                  q50 = quantile(diff.AI.mean, probs = 0.5, na.rm =T), 
                                                                  q75 = quantile(diff.AI.mean, probs = 0.75, na.rm =T), 
                                                                  q90 = quantile(diff.AI.mean, probs = 0.9, na.rm =T)) 

ano.quantiles <- cmip6s %>% summarise(q10 = quantile(diff.AI.mean, probs = 0.1, na.rm =T), 
                                      q25 = quantile(diff.AI.mean, probs = 0.25, na.rm =T), 
                                      q50 = quantile(diff.AI.mean, probs = 0.5, na.rm =T), 
                                      q75 = quantile(diff.AI.mean, probs = 0.75, na.rm =T), 
                                      q90 = quantile(diff.AI.mean, probs = 0.9, na.rm =T)) 
qbreaks <- c(-5, -1,-0.1, -0.01, 0, 0.01, 0.1, 1,5)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")

“diff.AI.mean” is the multimodel AI average minus the reference AI (period 1970-2000). More blue: AI mean is superior to AI ref, hence wetter than the ref. More red: the anomalie is negative, hence AI mean is inferior to AI ref, hence dryer.

AS OF NOW THE SAME QUANTILES ARE USED FOR ALL MAPS.

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  #quant = subset(ano.quantiles, model == "historical" & period == i)
  #b = quant[1,c(3:7)] %>% as.numeric
  b = qbreaks
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = diff.AI.mean))+
    borders(colour = "grey60")+
           geom_point(data = subset(cmip6s, period == i & model == "historical" & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s,  period == i & model == "historical"  & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
    binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale),
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    b = qbreaks
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = diff.AI.mean))+
      borders(colour = "grey60")+
                 geom_point(data = subset(cmip6s, period == i & model == "historical" & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s,  period == i & model == "historical"  & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
      binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale),
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Compared to 1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}



ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

6.5 Wetter or dryer, binned, % (reference period: 1970-2000)

cmip6s$diff.AI.mean.percent <- with(cmip6s, (AI.mean-AI.ref.mean)/AI.ref.mean*100)


pbreaks <- c(-40, -30,-20, -10, 0, 10, 20, 30, 40)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef","white","white", "#fddbc7", "#f4a582", "#d6604d", "#b2182b")

“diff.AI.mean” is the multimodel AI average minus the reference AI (period 1970-2000). More blue: AI mean is superior to AI ref, hence wetter than the ref. More red: the anomalie is negative, hence AI mean is inferior to AI ref, hence dryer.

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  #quant = subset(ano.quantiles, model == "historical" & period == i)
  #b = quant[1,c(3:7)] %>% as.numeric
  b = pbreaks
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = diff.AI.mean.percent))+
    borders(colour = "grey60")+
           geom_point(data = subset(cmip6s, period == i & model == "historical" & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s,  period == i & model == "historical"  & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
    binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale))+
    labs(title = i, fill = "")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom") %>% annotate_figure(top = "% change of AI index compared to 1970-2000")

#,guide = guide_legend(label.theme = element_text(angle = 0))

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    b = pbreaks
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = diff.AI.mean.percent))+
      borders(colour = "grey60", size = 0.5)+
                 geom_point(data = subset(cmip6s, period == j & model == i & plus.AI > 10), aes(x=lon, y = lat), shape = "+", col = "grey20", size = 0.5)+
       geom_point(data = subset(cmip6s,  period == j & model == i & minus.AI > 10), aes(x=lon, y = lat), shape = "-", col = "grey20", size = 0.5)+
      binned_scale(aesthetics = "fill", breaks = b, palette = function(x) rev(colscale))+
      labs(fill = "% change of AI index\ncompared to 1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}


ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

7 Difference between future and reference category of AI, using cat maj

7.1 Detailed changes and colors cat maj

cmip6s$change.catmaj <- with(cmip6s, ifelse(cat.maj == cat.maj.ref, 0, 1))

cmip6s$change.cat.maj <- with(cmip6s, paste(cat.maj.ref, cat.maj, sep = " to "))

cmip6s$change.cat.maj <- with(cmip6s, ifelse(cat.maj.ref == "Arid" & cat.maj %in% c("Cold","Humid"),"Arid to cold or humid",
         ifelse(cat.maj.ref == "Arid" & cat.maj %in% c("Semi-Arid","Dry subhumid"), "Arid to semi-arid or dry subhumid",
         ifelse(cat.maj.ref == "Arid" & cat.maj == "Hyperarid", "Arid to hyperarid",
         ifelse(cat.maj.ref == "Semi-arid" & cat.maj == "Dry subhumid", "Semi-arid to dry subhumid",
         ifelse(cat.maj.ref == "Semi-arid" & cat.maj == "Arid", "Semi-arid to arid",
         ifelse(cat.maj.ref == "Semi-arid" & cat.maj == "Hyperarid", "Semi-arid to hyperarid",
         ifelse(cat.maj.ref == "Semi-arid" & cat.maj %in% c("Humid","Cold"), "Semi-arid to humid or cold",
         ifelse(cat.maj.ref == "Dry subhumid" & cat.maj %in% c("Semi-arid","Arid","Hyperarid"), "Dry subhumid to arid, semi-arid or hyperarid",
         ifelse(cat.maj.ref == "Dry subhumid" & cat.maj %in% c("Humid","Cold"), "Dry subhumid to humid or cold",
         ifelse(cat.maj.ref == "Humid" & cat.maj == "Cold","Humid to cold",
         ifelse(cat.maj.ref == "Humid" & cat.maj %in% c("Dry subhumid", "Semi-arid"), "Humid to dry-subhumid or semi-arid",
         ifelse(cat.maj.ref == "Humid" & cat.maj %in% c("Arid", "Hyperarid"), "Humid to arid or hyperarid",
         ifelse(cat.maj.ref == "Cold" & cat.maj == "Humid", "Cold to humid", 
         ifelse(cat.maj.ref == "Cold" & cat.maj %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "Cold to dryland",
         ifelse(cat.maj.ref == "Hyperarid" & cat.maj == "Arid", "Hyperarid to arid",
         ifelse(cat.maj.ref == "Hyperarid" & cat.maj %in% c("Semi-arid","Dry-subhumid"),"Hyperarid to semi arid or dry subhumid",
         ifelse(cat.maj.ref == "Hyperarid" & cat.maj %in% c("Humid", "Cold"), "Hyperarid to humid or cold", NA))))))))))))))))))

cmip6s$change.cat.maj <- factor(cmip6s$change.cat.maj, 
                                levels = c("Arid to hyperarid","Semi-arid to arid", "Dry subhumid to arid, semi-arid or hyperarid", "Humid to dry-subhumid or semi-arid", "Cold to humid","Hyperarid to arid", "Arid to semi-arid or dry subhumid", "Dry subhumid to humid or cold", "Semi-arid to dry subhumid", "Semi-arid to humid or cold"))
table(cmip6s$change.cat.maj)
## 
##                            Arid to hyperarid 
##                                            0 
##                            Semi-arid to arid 
##                                            0 
## Dry subhumid to arid, semi-arid or hyperarid 
##                                            0 
##           Humid to dry-subhumid or semi-arid 
##                                            0 
##                                Cold to humid 
##                                            0 
##                            Hyperarid to arid 
##                                            0 
##            Arid to semi-arid or dry subhumid 
##                                            0 
##                Dry subhumid to humid or cold 
##                                            0 
##                    Semi-arid to dry subhumid 
##                                            0 
##                   Semi-arid to humid or cold 
##                                            0

colors_changes <- c("Hyperarid to arid" = "#E7D7B6", 
                    "Arid to semi-arid or dry subhumid" = colorvec[7],
                    "Dry subhumid to humid or cold" = colorvec[8],
                    "Semi-arid to dry subhumid" = colorvec[9],
                    "Semi-arid to humid or cold" = colorvec[10],
                    "Humid to cold" = colorvec[11],
                    "Cold to humid" = "#E7D7B6", 
                    "Humid to dry-subhumid or semi-arid" = colorvec[5],
                    "Dry subhumid to arid, semi-arid or hyperarid" = colorvec[4],
                    "Semi-arid to arid" = colorvec[3],
                    "Arid to hyperarid" = colorvec[1])

7.2 Changes towards dryer categories cat maj

7.2.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + 
    geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & cat.maj.dryer == "dryer"), 
                aes(x=lon, y = lat,  fill = change.cat.maj))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
  map_list[[i]] <- g
}


ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

7.2.1.1 Future

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmaj == 1 & cat.maj.dryer == "dryer"), 
                              aes(x=lon, y = lat,  fill = change.cat.maj))+
      #   geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & nb.maj > 10), aes(x= lon, y = lat, fill = "'"))+
      borders(colour = "grey60")+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      #  scale_fill_viridis_d()+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

7.3 Changes towards wetter categories cat maj

7.3.0.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & cat.maj.dryer == "wetter"), 
                              aes(x=lon, y = lat,  fill = change.cat.maj))+
    borders(colour = "grey60", ylim = c(-40,40))+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+
    ylim(-50,50)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

7.3.0.2 Future

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmaj == 1 & cat.maj.dryer == "wetter"), 
                              aes(x=lon, y = lat,  fill = change.cat.maj))+
      borders(colour = "grey60", ylim = c(-40,40))+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+ylim(-50,50)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list, nrow = 3, ncol = 2, common.legend = T, legend = "bottom")

8 Difference between future and reference category of AI, using cat mean

8.1 Detailed changes and colors

cmip6s$change.catmean <- with(cmip6s, ifelse(cat.AI == cat.AI.ref, 0, 1))

cmip6s$change.cat.mean <- with(cmip6s, paste(cat.AI.ref, cat.AI, sep = " to "))
unique(cmip6s$change.cat)
## NULL

# cmip6s$change.cat.mean <- with(cmip6s, 
#       ifelse(cat.AI.ref == "Arid" & cat.AI %in% c("Cold","Humid"),"Arid to cold or humid",
#       ifelse(cat.AI.ref == "Arid" & cat.AI %in% c("Semi-Arid","Dry subhumid"), "Arid to semi-arid or dry subhumid",
#       ifelse(cat.AI.ref == "Arid" & cat.AI == "Hyperarid", "Arid to hyperarid",
#       ifelse(cat.AI.ref == "Semi-arid" & cat.AI == "Dry subhumid", "Semi-arid to dry subhumid",
#       ifelse(cat.AI.ref == "Semi-arid" & cat.AI == "Arid", "Semi-arid to arid",
#       ifelse(cat.AI.ref == "Semi-arid" & cat.AI == "Hyperarid", "Semi-arid to hyperarid",
#       ifelse(cat.AI.ref == "Semi-arid" & cat.AI %in% c("Humid","Cold"), "Semi-arid to humid or cold",
#       ifelse(cat.AI.ref == "Dry subhumid" & cat.AI %in% c("Semi-arid","Arid","Hyperarid"), "Dry subhumid to arid, semi-arid or hyperarid",
#       ifelse(cat.AI.ref == "Dry subhumid" & cat.AI %in% c("Humid","Cold"), "Dry subhumid to humid or cold",
#       ifelse(cat.AI.ref == "Humid" & cat.AI == "Cold","Humid to cold",
#       ifelse(cat.AI.ref == "Humid" & cat.AI %in% c("Dry subhumid","Semi-arid"), "Humid to dry-subhumid or semi-arid",
#       ifelse(cat.AI.ref == "Humid" & cat.AI %in% c("Arid", "Hyperarid"), "Humid to arid or hyperarid",
#       ifelse(cat.AI.ref == "Cold" & cat.AI == "Humid", "Cold to humid", 
#       ifelse(cat.AI.ref == "Cold" & cat.AI %in% c("Dry subhumid","Semi-arid","Arid","Hyperarid"), "Cold to dryland",
#       ifelse(cat.AI.ref == "Hyperarid" & cat.AI == "Arid", "Hyperarid to arid",
#       ifelse(cat.AI.ref == "Hyperarid" & cat.AI %in% c("Semi-arid","Dry-subhumid"),"Hyperarid to semi arid or dry subhumid",
#       ifelse(cat.AI.ref == "Hyperarid" & cat.AI %in% c("Humid", "Cold"), "Hyperarid to humid or cold", 
#       NA))))))))))))))))))

table(cmip6s$change.cat.mean)
## 
##                 Arid to Arid            Arid to Hyperarid 
##                        11080                          639 
##            Arid to Semi-arid                 Cold to Cold 
##                          340                        95966 
## Dry subhumid to Dry subhumid        Dry subhumid to Humid 
##                         4022                          288 
##    Dry subhumid to Semi-arid        Humid to Dry subhumid 
##                         2242                         2484 
##               Humid to Humid           Humid to Semi-arid 
##                        57318                          198 
##            Hyperarid to Arid       Hyperarid to Hyperarid 
##                          240                         6402 
##                     NA to NA            Semi-arid to Arid 
##                         5130                         1209 
##    Semi-arid to Dry subhumid       Semi-arid to Semi-arid 
##                          277                        12514
cmip6s$change.cat.mean <- factor(cmip6s$change.cat.mean, 
                                  levels = c("Arid to Hyperarid","Semi-arid to Arid", "Dry subhumid to Semi-arid", "Humid to Semi-arid", "Humid to Dry subhumid", "Cold to Humid","Hyperarid to arid", "Arid to Semi-arid", "Dry subhumid to Humid", "Semi-arid to Dry subhumid", "Semi-arid to Humid"))

colors_changes <- c("Hyperarid to Arid" = "#E7D7B6", 
                     "Arid to Semi-arid" = colorvec[7],
                     "Dry subhumid to Humid" = colorvec[8],
                     "Semi-arid to Dry subhumid" = colorvec[9],
                     "Semi-arid to Humid" = colorvec[10],
                     "Humid to Cold" = colorvec[11],
                     "Cold to Humid" = "#E7D7B6",
                     "Humid to Dry subhumid" = colorvec[5],
                     "Humid to Semi-arid" = colorvec[4],
                     "Dry subhumid to Semi-arid" = colorvec[3],
                     "Semi-arid to Arid" = colorvec[2],
                     "Arid to Hyperarid" = colorvec[1])

8.2 Changes towards dryer categories cat mean

8.2.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + 
    geom_raster(data = subset(cmip6s, period == i & model == "historical" & change.catmean == 1 & cat.dryer == "dryer"), 
                aes(x=lon, y = lat,  fill = change.cat.mean))+
    borders(colour = "grey60")+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

8.2.1.1 Future

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmean == 1 & cat.dryer == "dryer"), 
                              aes(x=lon, y = lat,  fill = change.cat.mean))+
      borders(colour = "grey60")+
      scale_fill_manual(values = colors_changes)+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggarrange(plotlist = map_list, ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

8.3 Changes towards wetter categories cat mean

8.3.1 Historical

map_list <- list()

for(i in c("1850_1880","1985_2015")){
  g <- ggplot() + geom_tile(data = subset(cmip6s, period == i & model == "historical" & change.catmaj == 1 & cat.dryer == "wetter"), 
                            aes(x=lon, y = lat,  fill = change.cat.mean))+
    borders(colour = "grey60", ylim = c(-40,40))+
    scale_fill_manual(values =  colors_changes, na.translate = F)+
    labs(title = i, fill = "Change of category\ncompared to\n1970-2000")+
    theme_void()+ylim(-50,50)+
    theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
  map_list[[i]] <- g
}

ggarrange(plotlist = map_list, ncol = 2, common.legend = T, legend = "bottom")

8.3.2 Future

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i & change.catmean == 1 & cat.dryer == "wetter"), 
                              aes(x=lon, y = lat,  fill = change.cat.mean))+
      borders(colour = "grey60", ylim = c(-40,40))+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      labs(fill = "Change of category\ncompared to\n1970-2000", title = index)+
      theme_void()+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

8.4 Big map of SSP585 2070-2100

g <- ggplot() + 
  geom_tile(data = subset(cmip6s, period == "2070_2100" & model == "SSP585" & change.catmean == 1), aes(x=lon, y = lat,  fill = change.cat.mean))+
      borders(colour = "grey60")+
      scale_fill_manual(values = colors_changes, na.translate = F)+
      labs(fill = "Change of aridity category\nfor the period 2070-2100\nand SSP 5-8.5,\ncompared to 1970-2000", title = "")+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")+guides(fill = guide_legend(nrow = 3))
   
print(g)

9 Tables

9.1 Changes of aridity category using cat mean

tab <- cmip6s %>% group_by(model, period, Continent) %>% summarise(Cold = table(cat.AI)[1], Hyperarid = table(cat.AI)[2], Arid = table(cat.AI)[3], Semiarid = table(cat.AI)[4], Drysubhumid = table(cat.AI)[5], Humid = table(cat.AI)[6], total = n()) 
write.table(tab, "tab.catAI.txt")

tab.future <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(period, model, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(period + model ~cat.AI) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

kable(tab.future[order(tab.future$model, decreasing = F),]) %>% kable_styling(bootstrap_options = "bordered") %>%
  column_spec(c(7,10), italic = T, include_thead = T) %>% row_spec(2, bold = T) 

	model	period	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
1	historical	1850_1880	3.4	6.0	7.1	3.3	19.8	28.8	48.9	2.6
2	historical	1970_2000	3.3	6.0	7.0	3.3	19.6	27.8	49.9	2.6
3	historical	1985_2015	3.3	6.1	6.9	3.3	19.6	27.3	50.6	2.6
4	SSP245	2030_2060	3.5	6.3	7.8	3.3	20.9	28.8	47.7	2.6
7	SSP245	2070_2100	3.6	6.4	7.8	3.5	21.3	29.4	46.6	2.6
5	SSP370	2030_2060	3.2	6.7	7.6	3.0	20.5	28.8	48.1	2.6
8	SSP370	2070_2100	4.0	6.3	8.3	3.7	22.3	28.5	46.6	2.6
6	SSP585	2030_2060	3.6	6.2	7.9	3.5	21.2	29.3	47.0	2.6
9	SSP585	2070_2100	3.8	6.4	8.4	3.6	22.2	30.1	45.1	2.6

# ggplot(tab.future, aes(x = period))+
#   geom_point(aes(y = Hyperarid, col = "Hyperarid", shape = model))+
#   geom_line(aes(y=Hyperarid, col = "Hyperarid", group = model, lty = model))+
#   geom_point(aes(y = Arid, col = "Arid", shape = model))+
#   geom_line(aes(y=Arid, col = "Arid", group = model, lty = model))+
#   geom_point(aes(y = get('Semi-arid'), col = "Semi-arid", shape = model))+
#   geom_line(aes(y=get('Semi-arid'), col = "Semi-arid", group = model, lty = model))+
#   geom_point(aes(y = get('Dry subhumid'), col = "Dry subhumid", shape = model))+
#   geom_line(aes(y=get('Dry subhumid'), col = "Dry subhumid", group = model, lty = model))+
#   geom_point(aes(y = Humid, col = "Humid", shape = model))+
#   geom_line(aes(y= Humid, col = "Humid", group = model, lty = model))+
#   scale_color_manual(values = col.cat)+
#   theme_minimal()

cmip6$cat.AI <- factor(cmip6$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))
tab.percent <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() %>% 
  group_by(period, model, cat.AI) %>% 
  summarise(mmmean = mean(count, na.rm = T), mmsd = sd(count, na.rm = T)) %>%
  mutate(percent = round(mmmean/sum(mmmean)*100, 1), sd.percent = round(mmsd/sum(mmmean)*100,1)) 
write.table(tab.percent, "tab.percent.txt")

tab.percent <- read.table("tab.percent.txt")

df245 <- rbind(subset(tab.percent, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.percent, model == "SSP245")) %>%
  group_by(period, model) %>% 
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

df245$cat.AI <- factor(df245$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

b245 <- ggplot(data = df245, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = paste(percent, " ± ", sd.percent, sep = "")), size = 3)+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 2-4.5", y = "%", x = "", col = "", fill = "")+
  theme_minimal()

df370 <- rbind(subset(tab.percent, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.percent, model == "SSP370")) %>%
  group_by(period, model) %>% 
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

df370$cat.AI <- factor(df370$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

b370 <- ggplot(data = df370, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = paste(percent, " ± ", sd.percent, sep = "")), size = 3)+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 3-7.0", y = "%", x = "", col = "", fill = "")+
  theme_minimal()


df585 <- rbind(subset(tab.percent, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.percent, model == "SSP585")) %>%
  group_by(period, model) %>% 
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

df585$cat.AI <- factor(df585$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

b585 <- ggplot(data = df585, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = paste(percent, " ± ", sd.percent, sep = "")), size = 3)+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 5-8.5", y = "%", x = "", col = "", fill = "")+
  theme_minimal()


ggarrange(plotlist = list(b245, b370, b585), common.legend = T, legend = "bottom", ncol = 3)

9.2 Changes of aridity category using cat maj

tab.maj <- cmip6s %>% group_by(model, period, Continent) %>% summarise(Cold = table(cat.maj)[1], Hyperarid = table(cat.maj)[2], Arid = table(cat.maj)[3], Semiarid = table(cat.maj)[4], Drysubhumid = table(cat.maj)[5], Humid = table(cat.maj)[6], total = n()) 
write.table(tab.maj, "tab.catmaj.txt")

tab.future.maj <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC")) %>%
  group_by(period, model, cat.maj) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL) %>%
  reshape2::dcast(period + model ~cat.maj) %>%
  mutate("Sum drylands" = rowSums(.[c("Hyperarid","Arid","Semi-arid","Dry subhumid")])) %>%
  select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))

kable(tab.future.maj[order(tab.future.maj$model, decreasing = F),]) %>% kable_styling(bootstrap_options = "bordered") %>%
  column_spec(c(7,10), italic = T, include_thead = T) %>% row_spec(2, bold = T) 

	model	period	Hyperarid	Arid	Semi-arid	Dry subhumid	Sum drylands	Humid	Cold	NA
1	historical	1850_1880	4.2	7.1	8.7	0.8	20.8	28.3	48.3	2.6
2	historical	1970_2000	4.1	7.1	8.6	0.7	20.5	28.1	48.9	2.6
3	historical	1985_2015	4.0	7.0	8.4	0.9	20.3	27.3	49.9	2.6
4	SSP245	2030_2060	4.0	7.8	9.3	0.8	21.9	28.4	47.2	2.6
7	SSP245	2070_2100	4.2	7.7	9.4	0.9	22.2	29.1	46.1	2.6
5	SSP370	2030_2060	4.1	7.2	9.4	0.8	21.5	28.6	47.3	2.6
8	SSP370	2070_2100	4.5	7.7	10.0	0.9	23.1	28.4	45.9	2.6
6	SSP585	2030_2060	4.2	7.5	9.6	1.0	22.3	28.8	46.3	2.6
9	SSP585	2070_2100	4.5	7.9	9.4	0.9	22.7	30.0	44.6	2.6

tab.flow.catmaj <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, cat.maj) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL)

df245 <- rbind(subset(tab.flow.catmaj, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow.catmaj, model == "SSP245")) %>%
  mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b245 <- ggplot(data = df245, aes(x = period, y = percent))+
  geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 2-4.5", y = "%", x = "")+
  theme_minimal()

df370 <- rbind(subset(tab.flow.catmaj, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow.catmaj, model == "SSP370")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b370 <- ggplot(data = df370, aes(x = period, y = percent))+
  geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 3-7.0", y = "%", x = "")+
  theme_minimal()


df585 <- rbind(subset(tab.flow.catmaj, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow.catmaj, model == "SSP585")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.maj %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b585 <- ggplot(data = df585, aes(x = period, y = percent))+
  geom_col(aes(group = cat.maj, col = cat.maj, fill = cat.maj))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  labs(title = "SSP 5-8.5", y = "%", x = "")+
  theme_minimal()


ggarrange(plotlist = list(b245, b370, b585), common.legend = T, legend = "bottom", ncol = 3)

9.3 Changes by continent, cat mean

cmip6$cat.AI <- factor(cmip6$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

tab.percent.cont <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, Continent, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() %>% 
  group_by(period, model, Continent, cat.AI) %>% 
  summarise(mmmean = mean(count, na.rm = T), mmsd = sd(count, na.rm = T)) %>%
  mutate(percent = round(mmmean/sum(mmmean)*100, 1), sd.percent = round(mmsd/sum(mmmean)*100,1)) 
write.table(tab.percent.cont, "tab.percent.continent.txt")

tab.percent.cont <- read.table("tab.percent.continent.txt")

tab_list_percent <- list()
tab_list_sd <- list()

for(i in unique(tab.percent.cont$Continent)){
  tab.i <- subset(tab.percent.cont, Continent == i)
  
  df.percent <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "percent")
  
  df.sd <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "sd.percent")
  # sd for a sum: square root of sum of squared sd
  
  drycats <- which(names(df.percent) %in% c("Hyperarid","Arid","Semi-arid","Dry subhumid"))
  
  df.percent <- df.percent %>% mutate("Sum drylands" = rowSums(.[drycats], na.rm = T))
  df.sd <- df.sd %>% mutate("Sum drylands" = round(sqrt(rowSums((.[drycats])^2, na.rm = T)),2))
  
  missing <- setdiff(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"), names(df.percent))
  
  df.percent[, missing] <- 0
  df.sd[,missing] <- 0
  
  df.percent <- df.percent %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.percent[is.na(df.percent)] <- 0
  
  tab_list_percent[[i]] <- df.percent
  
  df.sd <- df.sd %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.sd[is.na(df.sd)] <- 0
  
  tab_list_sd[[i]] <- df.sd
}

k_list <- list()


kdiff_list <- list()

for(i in names(tab_list_percent)){
  
  df <- tab_list_percent[[i]] %>% select(c("model", "period"))
  
  
  
  prop.ref <- tab_list_percent[[i]] %>% filter(period == "1970_2000")
  prop.ref <- prop.ref %>% rbind(prop.ref[rep(1,dim(df)[1]-1),])
  df2 <- tab_list_percent[[i]][,c(3:10)] - prop.ref[,c(3:10)] %>% round(1)
  diff.prop.percent <- cbind(df, df2)
  
  sd.ref <- tab_list_sd[[i]] %>% filter(period == "1970_2000")
  sd.ref <- sd.ref %>% rbind(sd.ref[rep(1,dim(df)[1]-1),])
  df3 <- sqrt(tab_list_sd[[i]][,c(3:10)]^2 + sd.ref[,c(3:10)]^2) %>% round(1)
  diff.sd.percent <- cbind(df, df3)
  
  df$Hyperarid <- paste(tab_list_percent[[i]]$Hyperarid, tab_list_sd[[i]]$Hyperarid, sep = " ± ")
  df$Arid <- paste(tab_list_percent[[i]]$Arid, tab_list_sd[[i]]$Arid, sep = " ± ")
  df$'Semi-Arid' <- paste(tab_list_percent[[i]][,5], tab_list_sd[[i]][,5], sep = " ± ")
  df$'Dry subhumid' <- paste(tab_list_percent[[i]][,6], tab_list_sd[[i]][,6], sep = " ± ")
  df$'Sum drylands' <- paste(tab_list_percent[[i]][,7], tab_list_sd[[i]][,7], sep = " ± ")
  df$Humid <- paste(tab_list_percent[[i]]$Humid, tab_list_sd[[i]]$Humid, sep = " ± ")
  df$Cold <- paste(tab_list_percent[[i]]$Cold, tab_list_sd[[i]]$Cold, sep = " ± ")
  
  df.diff <- df
  df.diff$Hyperarid <- paste(round(diff.prop.percent$Hyperarid, 1), round(diff.sd.percent$Hyperarid, 1), sep = " ± ")
  df.diff$Arid <- paste(round(diff.prop.percent$Arid, 1), round(diff.sd.percent$Arid, 1), sep = " ± ")
  df.diff$'Semi-Arid' <- paste(round(diff.prop.percent[,5], 1), round(diff.sd.percent[,5], 1), sep = " ± ")
  df.diff$'Dry subhumid' <- paste(round(diff.prop.percent[,6], 1), round(diff.sd.percent[,6], 1), sep = " ± ")
  df.diff$'Sum drylands' <- paste(round(diff.prop.percent[,7],1), round(diff.sd.percent[,7], 1), sep = " ± ")
  df.diff$Humid <- paste(round(diff.prop.percent$Humid, 1), round(diff.sd.percent$Humid, 1), sep = " ± ")
  df.diff$Cold <- paste(round(diff.prop.percent$Cold, 1), round(diff.sd.percent$Cold, 1), sep = " ± ")
  
  k_list[[i]] <- kable(df[order(df$model, decreasing = F),], caption = i) %>% kable_styling(bootstrap_options = "bordered") %>%
    column_spec(8, italic = T, background = colorvec[5]) %>% row_spec(2, bold = T) 
  kdiff_list[[i]] <- kable(df.diff[order(df.diff$model, decreasing = F),], caption = i) %>% kable_styling(bootstrap_options = "bordered") %>% 
    column_spec(8, italic = T, background = colorvec[5]) %>% row_spec(2, bold = T) 
  
}

9.3.1 Proportion of aridity categories by continent

k_list

$AFRICA

AFRICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	27.2 ± 4.7	17.1 ± 4.9	15.6 ± 3.6	6.3 ± 1.4	66.2 ± 7.81	33.7 ± 7.1	0 ± 0
2	historical	1970_2000	26.9 ± 4.7	16.7 ± 4.2	15.1 ± 3.4	6.2 ± 1.5	64.9 ± 7.32	32.8 ± 6.2	2.3 ± 0
3	historical	1985_2015	26.5 ± 4.8	17.1 ± 4.6	15.4 ± 3.3	6 ± 1.5	65 ± 7.57	32.7 ± 6.3	2.2 ± 0
4	SSP245	2030_2060	27 ± 4.8	18 ± 4.9	16.5 ± 3.5	6.4 ± 1.7	67.9 ± 7.89	32 ± 7.6	0 ± 0
7	SSP245	2070_2100	27.9 ± 4.7	17.4 ± 5	17.2 ± 3.4	6.2 ± 1.5	68.7 ± 7.8	31.3 ± 7.9	0 ± 0
5	SSP370	2030_2060	26 ± 5	16.6 ± 3.6	17.2 ± 1.8	6.8 ± 1	66.6 ± 6.5	32.3 ± 6.7	1.1 ± 0
8	SSP370	2070_2100	28 ± 3.5	17.6 ± 4.8	16.6 ± 3.1	6.2 ± 1.6	68.4 ± 6.89	30.4 ± 6.4	1.2 ± 0
6	SSP585	2030_2060	27.1 ± 4.6	17.6 ± 4.8	16.8 ± 3.5	6.3 ± 1.5	67.8 ± 7.66	32.1 ± 7.5	0 ± 0
9	SSP585	2070_2100	29.7 ± 10.5	17.3 ± 5.3	16.2 ± 4.2	6 ± 1.5	69.2 ± 12.58	30.9 ± 8.6	0 ± 0

$ASIA

ASIA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	3.9 ± 2.3	12.2 ± 1.9	9.2 ± 1.5	4 ± 1.4	29.3 ± 3.62	32.4 ± 7.7	38 ± 9.2
2	historical	1970_2000	3.8 ± 2.1	12 ± 1.8	8.9 ± 1.5	3.8 ± 1.1	28.5 ± 3.33	30.7 ± 8.5	40.5 ± 9.6
3	historical	1985_2015	3.7 ± 2.2	11.8 ± 1.6	8.6 ± 1.5	3.7 ± 1.1	27.8 ± 3.3	30.1 ± 7.6	41.8 ± 8.8
4	SSP245	2030_2060	3.5 ± 1.8	12.5 ± 1.6	9.4 ± 1.6	4.3 ± 1.7	29.7 ± 3.35	34.3 ± 8.1	35.8 ± 9.2
7	SSP245	2070_2100	3.7 ± 2	12.5 ± 1.7	9.5 ± 1.7	4.2 ± 1.7	29.9 ± 3.56	36.2 ± 8.6	33.7 ± 9.3
5	SSP370	2030_2060	3.4 ± 2.3	12.4 ± 1.5	8.9 ± 1.7	3.7 ± 1.3	28.4 ± 3.48	34 ± 7.2	37.3 ± 9.5
8	SSP370	2070_2100	4.6 ± 1.9	12.5 ± 1.5	9.6 ± 1.5	3.8 ± 1.2	30.5 ± 3.09	35.2 ± 9.2	34.1 ± 9.7
6	SSP585	2030_2060	3.8 ± 2	12.5 ± 1.6	9.5 ± 1.8	4.2 ± 1.6	30 ± 3.52	35.4 ± 8.4	34.4 ± 9.3
9	SSP585	2070_2100	4.8 ± 3.2	13.1 ± 3.2	10.4 ± 3.4	4.2 ± 1.5	32.5 ± 5.86	36.5 ± 11.7	30.7 ± 9.9

$CENTRAL-AMERICA

CENTRAL-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	8.4 ± 3.6	30.1 ± 5.1	16.1 ± 8.7	54.6 ± 10.71	45 ± 6.1	0 ± 0
2	historical	1970_2000	0 ± 0	8.4 ± 4.1	29.7 ± 7.6	14.7 ± 7.8	52.8 ± 11.64	44.9 ± 12	2 ± 0
3	historical	1985_2015	0 ± 0	9 ± 4.2	30.7 ± 7	14.8 ± 7.4	54.5 ± 11.02	43.8 ± 11.5	1.4 ± 0
4	SSP245	2030_2060	0 ± 0	11.9 ± 6.2	36.3 ± 8	13.9 ± 6.2	62.1 ± 11.87	37.6 ± 8.6	0 ± 0
7	SSP245	2070_2100	0 ± 0	12.3 ± 5.5	38.1 ± 8.5	14.4 ± 7.9	64.8 ± 12.84	34.8 ± 8.6	0 ± 0
5	SSP370	2030_2060	0 ± 0	12.6 ± 6.3	31.4 ± 7.1	15 ± 7.9	59 ± 12.35	40.7 ± 11.8	0 ± 0
8	SSP370	2070_2100	4 ± 5.2	16 ± 5.5	35.5 ± 9.5	13.2 ± 5.2	68.7 ± 13.21	30.9 ± 10.6	0 ± 0
6	SSP585	2030_2060	0 ± 0	13 ± 7	37.6 ± 9	14.1 ± 7.1	64.7 ± 13.43	35 ± 9.4	0 ± 0
9	SSP585	2070_2100	2 ± 2.4	16.6 ± 8	41.2 ± 7.7	13.3 ± 7.7	73.1 ± 13.72	26.6 ± 8.8	0 ± 0

$EUROPE

EUROPE

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	1.9 ± 1.2	9.6 ± 4.1	4.5 ± 2.8	16 ± 5.11	49 ± 15.3	34.8 ± 20.2
2	historical	1970_2000	0 ± 0	1.5 ± 1	8.9 ± 3.4	3.9 ± 2.1	14.3 ± 4.12	44.9 ± 16.1	40.6 ± 20.4
3	historical	1985_2015	0 ± 0	1.3 ± 1	8.4 ± 2.9	3.8 ± 2	13.5 ± 3.66	46 ± 16.1	40.3 ± 19.9
4	SSP245	2030_2060	0 ± 0	1.9 ± 0.9	10.3 ± 3.6	5.1 ± 1.9	17.3 ± 4.17	52.7 ± 15.4	29.8 ± 19.7
7	SSP245	2070_2100	0 ± 0	2.2 ± 1.1	10.1 ± 4.8	5.1 ± 2.3	17.4 ± 5.44	55.4 ± 14.4	27 ± 18.6
5	SSP370	2030_2060	0 ± 0	2.1 ± 1	9.7 ± 4.7	4.8 ± 2.5	16.6 ± 5.42	52.7 ± 15.7	30.4 ± 20
8	SSP370	2070_2100	0 ± 0	3.1 ± 1.1	11.4 ± 5.9	6.4 ± 3.8	20.9 ± 7.1	52.8 ± 14.3	26 ± 18.6
6	SSP585	2030_2060	0 ± 0	2.3 ± 1.1	10.5 ± 3.7	4.9 ± 2.4	17.7 ± 4.55	55.2 ± 15.8	26.9 ± 19.9
9	SSP585	2070_2100	0 ± 0	3 ± 1.7	11.1 ± 5.4	6.1 ± 2.9	20.2 ± 6.36	57.2 ± 12.8	22.5 ± 17.3

$EUROPE-AFRICA

EUROPE-AFRICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	5.9 ± 4.9	23.6 ± 7	21.3 ± 8.8	8.8 ± 3.3	59.6 ± 12.7	24.9 ± 14.2	14.5 ± 0
2	historical	1970_2000	7 ± 5.5	26.5 ± 7.5	23.3 ± 10.3	9.9 ± 3.2	66.7 ± 14.24	27.5 ± 16.5	4.5 ± 7.5
3	historical	1985_2015	7.1 ± 5.9	26.3 ± 7.1	22.2 ± 7.7	10.3 ± 4.2	65.9 ± 12.73	27.8 ± 14.3	5.2 ± 7
4	SSP245	2030_2060	9.9 ± 6.6	28.9 ± 8.6	26.7 ± 8.9	9.4 ± 3.4	74.9 ± 14.43	23.1 ± 16.2	0.9 ± 0
7	SSP245	2070_2100	11.7 ± 5.7	28.6 ± 8.8	27.2 ± 8.1	9.5 ± 3.1	77 ± 13.61	21.9 ± 13.9	0 ± 0
5	SSP370	2030_2060	10.3 ± 6.3	28.2 ± 8.4	26 ± 8.7	9 ± 2.8	73.5 ± 13.92	24.8 ± 14	0.5 ± 0
8	SSP370	2070_2100	13.2 ± 6.2	27.3 ± 8.6	28.3 ± 9.7	8.7 ± 3.6	77.5 ± 14.81	21.2 ± 16.2	0 ± 0
6	SSP585	2030_2060	11.1 ± 6.4	28.8 ± 8.3	27.5 ± 8.8	8.7 ± 3	76.1 ± 14.01	22.6 ± 14.1	0.2 ± 0
9	SSP585	2070_2100	15.2 ± 5.7	29.1 ± 10.1	27.9 ± 8.8	8 ± 2.8	80.2 ± 14.82	18.6 ± 13.6	0 ± 0

$NORTH-AMERICA

NORTH-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	0.4 ± 0.2	5.8 ± 5.6	4.3 ± 2.6	10.5 ± 6.18	35.1 ± 9	54.2 ± 10.9
2	historical	1970_2000	0 ± 0	0.4 ± 0.2	4.3 ± 4.4	4.2 ± 2.4	8.9 ± 5.02	34.5 ± 9.8	56.5 ± 11.7
3	historical	1985_2015	0 ± 0	0.3 ± 0.2	4.6 ± 3.8	4.1 ± 2.2	9 ± 4.4	33.6 ± 10.3	57.2 ± 11.5
4	SSP245	2030_2060	0 ± 0	0.4 ± 0.3	6.5 ± 5.5	4.7 ± 2.6	11.6 ± 6.09	37.5 ± 8.7	50.8 ± 11.1
7	SSP245	2070_2100	0 ± 0	0.5 ± 0.3	6.4 ± 5.1	5.1 ± 2.9	12 ± 5.87	39.7 ± 9	48.1 ± 11
5	SSP370	2030_2060	0 ± 0	0.5 ± 0.4	5.9 ± 4.8	4.6 ± 2.7	11 ± 5.52	37 ± 9.7	51.8 ± 12.2
8	SSP370	2070_2100	0.1 ± 0	0.8 ± 0.7	7.6 ± 4	5.2 ± 2	13.7 ± 4.53	38.4 ± 10.6	47.7 ± 11.2
6	SSP585	2030_2060	0 ± 0	0.5 ± 0.4	6.6 ± 5.2	5.2 ± 2.8	12.3 ± 5.92	38.2 ± 9.1	49.3 ± 11.7
9	SSP585	2070_2100	0.2 ± 0	1 ± 1.4	7 ± 6.1	5.2 ± 2.5	13.4 ± 6.74	42.1 ± 9.3	44.5 ± 11.3

$OCEANIA

OCEANIA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	5.3 ± 0	32 ± 26	36.7 ± 17.1	6.6 ± 3.1	80.6 ± 31.27	18.9 ± 12.5	0.3 ± 0.2
2	historical	1970_2000	5.2 ± 0	31.3 ± 26.4	36.8 ± 17.4	6.8 ± 3.1	80.1 ± 31.77	19.3 ± 12.8	0.5 ± 0.4
3	historical	1985_2015	5.2 ± 0	33.4 ± 26.1	35.5 ± 17.9	6.8 ± 3.2	80.9 ± 31.81	18.4 ± 12.3	0.6 ± 0.3
4	SSP245	2030_2060	10.4 ± 0	35.8 ± 23.5	31.7 ± 15.1	5.7 ± 2.6	83.6 ± 28.05	15.9 ± 11.9	0.3 ± 0.3
7	SSP245	2070_2100	10.5 ± 0	37.9 ± 22.2	29.1 ± 13.1	6.2 ± 3.3	83.7 ± 25.99	16.2 ± 12	0.1 ± 0
5	SSP370	2030_2060	10.7 ± 0	34.8 ± 22.6	33.4 ± 15.3	5.8 ± 2.5	84.7 ± 27.41	14.9 ± 8.1	0.3 ± 0.2
8	SSP370	2070_2100	10.5 ± 10.2	38.1 ± 19.7	28.5 ± 13.8	5.7 ± 3.2	82.8 ± 26.32	16.8 ± 16.3	0.3 ± 0.3
6	SSP585	2030_2060	5.6 ± 3.6	35.6 ± 25.3	33.8 ± 16.7	6.7 ± 3.6	81.7 ± 30.74	17.7 ± 13.8	0.3 ± 0.3
9	SSP585	2070_2100	9.3 ± 6.4	39.4 ± 21.7	29.4 ± 12.9	5.5 ± 2.5	83.6 ± 26.16	16.2 ± 12.5	0 ± 0

$SOUTH-AMERICA

SOUTH-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0.1 ± 0.1	2.1 ± 1.7	15.2 ± 8.7	7.4 ± 2.6	24.8 ± 9.24	73.2 ± 11.7	1.9 ± 1
2	historical	1970_2000	0.1 ± 0.1	1.8 ± 1.3	13.9 ± 4.3	7.3 ± 2.3	23.1 ± 5.05	74.5 ± 7.3	2.3 ± 2.6
3	historical	1985_2015	0.2 ± 0.1	2 ± 1.2	14.5 ± 5.1	7.5 ± 2.4	24.2 ± 5.76	73.6 ± 7.7	2.1 ± 2.6
4	SSP245	2030_2060	0.2 ± 0.2	2.7 ± 2.2	17.4 ± 9.6	8.8 ± 3.6	29.1 ± 10.49	69.4 ± 13.1	1.3 ± 0.9
7	SSP245	2070_2100	0.2 ± 0.1	3.3 ± 3.5	18.3 ± 9.1	9.3 ± 4.6	31.1 ± 10.78	67.5 ± 14.7	1.3 ± 0.7
5	SSP370	2030_2060	0.2 ± 0.1	2.2 ± 1.7	16.3 ± 5.5	8.8 ± 3.7	27.5 ± 6.84	70.5 ± 8.9	1.9 ± 1.6
8	SSP370	2070_2100	0.2 ± 0.1	3.3 ± 2.5	18.6 ± 7	8.7 ± 3.8	30.8 ± 8.35	67.3 ± 10.8	1.8 ± 1.8
6	SSP585	2030_2060	0.2 ± 0.1	3.1 ± 2.5	18.6 ± 9	9.2 ± 4.4	31.1 ± 10.33	67.6 ± 13.1	1.2 ± 0.8
9	SSP585	2070_2100	0.6 ± 1	4.4 ± 3.8	19.6 ± 9.5	9.7 ± 3.9	34.3 ± 11	64.4 ± 12.6	1.1 ± 0.7

9.3.2 Changes in proportion of aridity categories by continent

kdiff_list

$AFRICA

AFRICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0.3 ± 6.6	0.4 ± 6.5	0.5 ± 5	0.1 ± 2.1	1.3 ± 10.7	0.9 ± 9.4	-2.3 ± 0
2	historical	1970_2000	0 ± 6.6	0 ± 5.9	0 ± 4.8	0 ± 2.1	0 ± 10.4	0 ± 8.8	0 ± 0
3	historical	1985_2015	-0.4 ± 6.7	0.4 ± 6.2	0.3 ± 4.7	-0.2 ± 2.1	0.1 ± 10.5	-0.1 ± 8.8	-0.1 ± 0
4	SSP245	2030_2060	0.1 ± 6.7	1.3 ± 6.5	1.4 ± 4.9	0.2 ± 2.3	3 ± 10.8	-0.8 ± 9.8	-2.3 ± 0
7	SSP245	2070_2100	1 ± 6.6	0.7 ± 6.5	2.1 ± 4.8	0 ± 2.1	3.8 ± 10.7	-1.5 ± 10	-2.3 ± 0
5	SSP370	2030_2060	-0.9 ± 6.9	-0.1 ± 5.5	2.1 ± 3.8	0.6 ± 1.8	1.7 ± 9.8	-0.5 ± 9.1	-1.2 ± 0
8	SSP370	2070_2100	1.1 ± 5.9	0.9 ± 6.4	1.5 ± 4.6	0 ± 2.2	3.5 ± 10.1	-2.4 ± 8.9	-1.1 ± 0
6	SSP585	2030_2060	0.2 ± 6.6	0.9 ± 6.4	1.7 ± 4.9	0.1 ± 2.1	2.9 ± 10.6	-0.7 ± 9.7	-2.3 ± 0
9	SSP585	2070_2100	2.8 ± 11.5	0.6 ± 6.8	1.1 ± 5.4	-0.2 ± 2.1	4.3 ± 14.6	-1.9 ± 10.6	-2.3 ± 0

$ASIA

ASIA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0.1 ± 3.1	0.2 ± 2.6	0.3 ± 2.1	0.2 ± 1.8	0.8 ± 4.9	1.7 ± 11.5	-2.5 ± 13.3
2	historical	1970_2000	0 ± 3	0 ± 2.5	0 ± 2.1	0 ± 1.6	0 ± 4.7	0 ± 12	0 ± 13.6
3	historical	1985_2015	-0.1 ± 3	-0.2 ± 2.4	-0.3 ± 2.1	-0.1 ± 1.6	-0.7 ± 4.7	-0.6 ± 11.4	1.3 ± 13
4	SSP245	2030_2060	-0.3 ± 2.8	0.5 ± 2.4	0.5 ± 2.2	0.5 ± 2	1.2 ± 4.7	3.6 ± 11.7	-4.7 ± 13.3
7	SSP245	2070_2100	-0.1 ± 2.9	0.5 ± 2.5	0.6 ± 2.3	0.4 ± 2	1.4 ± 4.9	5.5 ± 12.1	-6.8 ± 13.4
5	SSP370	2030_2060	-0.4 ± 3.1	0.4 ± 2.3	0 ± 2.3	-0.1 ± 1.7	-0.1 ± 4.8	3.3 ± 11.1	-3.2 ± 13.5
8	SSP370	2070_2100	0.8 ± 2.8	0.5 ± 2.3	0.7 ± 2.1	0 ± 1.6	2 ± 4.5	4.5 ± 12.5	-6.4 ± 13.6
6	SSP585	2030_2060	0 ± 2.9	0.5 ± 2.4	0.6 ± 2.3	0.4 ± 1.9	1.5 ± 4.8	4.7 ± 12	-6.1 ± 13.4
9	SSP585	2070_2100	1 ± 3.8	1.1 ± 3.7	1.5 ± 3.7	0.4 ± 1.9	4 ± 6.7	5.8 ± 14.5	-9.8 ± 13.8

$CENTRAL-AMERICA

CENTRAL-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	0 ± 5.5	0.4 ± 9.2	1.4 ± 11.7	1.8 ± 15.8	0.1 ± 13.5	-2 ± 0
2	historical	1970_2000	0 ± 0	0 ± 5.8	0 ± 10.7	0 ± 11	0 ± 16.5	0 ± 17	0 ± 0
3	historical	1985_2015	0 ± 0	0.6 ± 5.9	1 ± 10.3	0.1 ± 10.8	1.7 ± 16	-1.1 ± 16.6	-0.6 ± 0
4	SSP245	2030_2060	0 ± 0	3.5 ± 7.4	6.6 ± 11	-0.8 ± 10	9.3 ± 16.6	-7.3 ± 14.8	-2 ± 0
7	SSP245	2070_2100	0 ± 0	3.9 ± 6.9	8.4 ± 11.4	-0.3 ± 11.1	12 ± 17.3	-10.1 ± 14.8	-2 ± 0
5	SSP370	2030_2060	0 ± 0	4.2 ± 7.5	1.7 ± 10.4	0.3 ± 11.1	6.2 ± 17	-4.2 ± 16.8	-2 ± 0
8	SSP370	2070_2100	4 ± 5.2	7.6 ± 6.9	5.8 ± 12.2	-1.5 ± 9.4	15.9 ± 17.6	-14 ± 16	-2 ± 0
6	SSP585	2030_2060	0 ± 0	4.6 ± 8.1	7.9 ± 11.8	-0.6 ± 10.5	11.9 ± 17.8	-9.9 ± 15.2	-2 ± 0
9	SSP585	2070_2100	2 ± 2.4	8.2 ± 9	11.5 ± 10.8	-1.4 ± 11	20.3 ± 18	-18.3 ± 14.9	-2 ± 0

$EUROPE

EUROPE

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	0.4 ± 1.6	0.7 ± 5.3	0.6 ± 3.5	1.7 ± 6.6	4.1 ± 22.2	-5.8 ± 28.7
2	historical	1970_2000	0 ± 0	0 ± 1.4	0 ± 4.8	0 ± 3	0 ± 5.8	0 ± 22.8	0 ± 28.8
3	historical	1985_2015	0 ± 0	-0.2 ± 1.4	-0.5 ± 4.5	-0.1 ± 2.9	-0.8 ± 5.5	1.1 ± 22.8	-0.3 ± 28.5
4	SSP245	2030_2060	0 ± 0	0.4 ± 1.3	1.4 ± 5	1.2 ± 2.8	3 ± 5.9	7.8 ± 22.3	-10.8 ± 28.4
7	SSP245	2070_2100	0 ± 0	0.7 ± 1.5	1.2 ± 5.9	1.2 ± 3.1	3.1 ± 6.8	10.5 ± 21.6	-13.6 ± 27.6
5	SSP370	2030_2060	0 ± 0	0.6 ± 1.4	0.8 ± 5.8	0.9 ± 3.3	2.3 ± 6.8	7.8 ± 22.5	-10.2 ± 28.6
8	SSP370	2070_2100	0 ± 0	1.6 ± 1.5	2.5 ± 6.8	2.5 ± 4.3	6.6 ± 8.2	7.9 ± 21.5	-14.6 ± 27.6
6	SSP585	2030_2060	0 ± 0	0.8 ± 1.5	1.6 ± 5	1 ± 3.2	3.4 ± 6.1	10.3 ± 22.6	-13.7 ± 28.5
9	SSP585	2070_2100	0 ± 0	1.5 ± 2	2.2 ± 6.4	2.2 ± 3.6	5.9 ± 7.6	12.3 ± 20.6	-18.1 ± 26.7

$EUROPE-AFRICA

EUROPE-AFRICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	-1.1 ± 7.4	-2.9 ± 10.3	-2 ± 13.5	-1.1 ± 4.6	-7.1 ± 19.1	-2.6 ± 21.8	10 ± 7.5
2	historical	1970_2000	0 ± 7.8	0 ± 10.6	0 ± 14.6	0 ± 4.5	0 ± 20.1	0 ± 23.3	0 ± 10.6
3	historical	1985_2015	0.1 ± 8.1	-0.2 ± 10.3	-1.1 ± 12.9	0.4 ± 5.3	-0.8 ± 19.1	0.3 ± 21.8	0.7 ± 10.3
4	SSP245	2030_2060	2.9 ± 8.6	2.4 ± 11.4	3.4 ± 13.6	-0.5 ± 4.7	8.2 ± 20.3	-4.4 ± 23.1	-3.6 ± 7.5
7	SSP245	2070_2100	4.7 ± 7.9	2.1 ± 11.6	3.9 ± 13.1	-0.4 ± 4.5	10.3 ± 19.7	-5.6 ± 21.6	-4.5 ± 7.5
5	SSP370	2030_2060	3.3 ± 8.4	1.7 ± 11.3	2.7 ± 13.5	-0.9 ± 4.3	6.8 ± 19.9	-2.7 ± 21.6	-4 ± 7.5
8	SSP370	2070_2100	6.2 ± 8.3	0.8 ± 11.4	5 ± 14.1	-1.2 ± 4.8	10.8 ± 20.5	-6.3 ± 23.1	-4.5 ± 7.5
6	SSP585	2030_2060	4.1 ± 8.4	2.3 ± 11.2	4.2 ± 13.5	-1.2 ± 4.4	9.4 ± 20	-4.9 ± 21.7	-4.3 ± 7.5
9	SSP585	2070_2100	8.2 ± 7.9	2.6 ± 12.6	4.6 ± 13.5	-1.9 ± 4.3	13.5 ± 20.6	-8.9 ± 21.4	-4.5 ± 7.5

$NORTH-AMERICA

NORTH-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0	0 ± 0.3	1.5 ± 7.1	0.1 ± 3.5	1.6 ± 8	0.6 ± 13.3	-2.3 ± 16
2	historical	1970_2000	0 ± 0	0 ± 0.3	0 ± 6.2	0 ± 3.4	0 ± 7.1	0 ± 13.9	0 ± 16.5
3	historical	1985_2015	0 ± 0	-0.1 ± 0.3	0.3 ± 5.8	-0.1 ± 3.3	0.1 ± 6.7	-0.9 ± 14.2	0.7 ± 16.4
4	SSP245	2030_2060	0 ± 0	0 ± 0.4	2.2 ± 7	0.5 ± 3.5	2.7 ± 7.9	3 ± 13.1	-5.7 ± 16.1
7	SSP245	2070_2100	0 ± 0	0.1 ± 0.4	2.1 ± 6.7	0.9 ± 3.8	3.1 ± 7.7	5.2 ± 13.3	-8.4 ± 16.1
5	SSP370	2030_2060	0 ± 0	0.1 ± 0.4	1.6 ± 6.5	0.4 ± 3.6	2.1 ± 7.5	2.5 ± 13.8	-4.7 ± 16.9
8	SSP370	2070_2100	0.1 ± 0	0.4 ± 0.7	3.3 ± 5.9	1 ± 3.1	4.8 ± 6.8	3.9 ± 14.4	-8.8 ± 16.2
6	SSP585	2030_2060	0 ± 0	0.1 ± 0.4	2.3 ± 6.8	1 ± 3.7	3.4 ± 7.8	3.7 ± 13.4	-7.2 ± 16.5
9	SSP585	2070_2100	0.2 ± 0	0.6 ± 1.4	2.7 ± 7.5	1 ± 3.5	4.5 ± 8.4	7.6 ± 13.5	-12 ± 16.3

$OCEANIA

OCEANIA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0.1 ± 0	0.7 ± 37.1	-0.1 ± 24.4	-0.2 ± 4.4	0.5 ± 44.6	-0.4 ± 17.9	-0.2 ± 0.4
2	historical	1970_2000	0 ± 0	0 ± 37.3	0 ± 24.6	0 ± 4.4	0 ± 44.9	0 ± 18.1	0 ± 0.6
3	historical	1985_2015	0 ± 0	2.1 ± 37.1	-1.3 ± 25	0 ± 4.5	0.8 ± 45	-0.9 ± 17.8	0.1 ± 0.5
4	SSP245	2030_2060	5.2 ± 0	4.5 ± 35.3	-5.1 ± 23	-1.1 ± 4	3.5 ± 42.4	-3.4 ± 17.5	-0.2 ± 0.5
7	SSP245	2070_2100	5.3 ± 0	6.6 ± 34.5	-7.7 ± 21.8	-0.6 ± 4.5	3.6 ± 41	-3.1 ± 17.5	-0.4 ± 0.4
5	SSP370	2030_2060	5.5 ± 0	3.5 ± 34.8	-3.4 ± 23.2	-1 ± 4	4.6 ± 42	-4.4 ± 15.1	-0.2 ± 0.4
8	SSP370	2070_2100	5.3 ± 10.2	6.8 ± 32.9	-8.3 ± 22.2	-1.1 ± 4.5	2.7 ± 41.3	-2.5 ± 20.7	-0.2 ± 0.5
6	SSP585	2030_2060	0.4 ± 3.6	4.3 ± 36.6	-3 ± 24.1	-0.1 ± 4.8	1.6 ± 44.2	-1.6 ± 18.8	-0.2 ± 0.5
9	SSP585	2070_2100	4.1 ± 6.4	8.1 ± 34.2	-7.4 ± 21.7	-1.3 ± 4	3.5 ± 41.2	-3.1 ± 17.9	-0.5 ± 0.4

$SOUTH-AMERICA

SOUTH-AMERICA

	model	period	Hyperarid	Arid	Semi-Arid	Dry subhumid	Sum drylands	Humid	Cold
1	historical	1850_1880	0 ± 0.1	0.3 ± 2.1	1.3 ± 9.7	0.1 ± 3.5	1.7 ± 10.5	-1.3 ± 13.8	-0.4 ± 2.8
2	historical	1970_2000	0 ± 0.1	0 ± 1.8	0 ± 6.1	0 ± 3.3	0 ± 7.1	0 ± 10.3	0 ± 3.7
3	historical	1985_2015	0.1 ± 0.1	0.2 ± 1.8	0.6 ± 6.7	0.2 ± 3.3	1.1 ± 7.7	-0.9 ± 10.6	-0.2 ± 3.7
4	SSP245	2030_2060	0.1 ± 0.2	0.9 ± 2.6	3.5 ± 10.5	1.5 ± 4.3	6 ± 11.6	-5.1 ± 15	-1 ± 2.8
7	SSP245	2070_2100	0.1 ± 0.1	1.5 ± 3.7	4.4 ± 10.1	2 ± 5.1	8 ± 11.9	-7 ± 16.4	-1 ± 2.7
5	SSP370	2030_2060	0.1 ± 0.1	0.4 ± 2.1	2.4 ± 7	1.5 ± 4.4	4.4 ± 8.5	-4 ± 11.5	-0.4 ± 3.1
8	SSP370	2070_2100	0.1 ± 0.1	1.5 ± 2.8	4.7 ± 8.2	1.4 ± 4.4	7.7 ± 9.8	-7.2 ± 13	-0.5 ± 3.2
6	SSP585	2030_2060	0.1 ± 0.1	1.3 ± 2.8	4.7 ± 10	1.9 ± 5	8 ± 11.5	-6.9 ± 15	-1.1 ± 2.7
9	SSP585	2070_2100	0.5 ± 1	2.6 ± 4	5.7 ± 10.4	2.4 ± 4.5	11.2 ± 12.1	-10.1 ± 14.6	-1.2 ± 2.7

9.3.3 Changes in proportion of aridity categories, by continent, calculated by model

cmip6$cat.AI <- factor(cmip6$cat.AI, levels = c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid", "Cold"))

tab.percent.cont <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, Continent, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() %>% 
  group_by(period, model, Continent, cat.AI) %>% 
  summarise(mmmean = mean(count, na.rm = T), mmsd = sd(count, na.rm = T)) %>%
  mutate(percent = round(mmmean/sum(mmmean)*100, 1), sd.percent = round(mmsd/sum(mmmean)*100,1)) 
write.table(tab.percent.cont, "tab.percent.continent.txt")

tab.cont <- cmip6 %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, Continent, cat.AI, source) %>%
  summarise(count = n()) %>%
  ungroup() 

list_ref <- list()

for(i in unique(tab.cont$source)){
  tab.cont.i <- filter(tab.cont, source == i)
  tab.cont.i.ref <- filter(tab.cont.i, period == "1970_2000")
  
  for(j in unique(tab.cont.i$Continent)){
    tab.cont.i.j <- filter(tab.cont.i, Continent == j)
    tab.cont.i.j.ref <- filter(tab.cont.i.ref, Continent == j)
    
    for(k in unique(tab.cont.i$cat.AI)){
      
      index <- paste(i, j, k, sep = "_")
      tab.cont.i.j.k <- filter(tab.cont.i.j, cat.AI == k)
      
      lines <- ifelse(dim(tab.cont.i.j.k)[1] > 0, dim(tab.cont.i.j.k)[1] - 1, 1)
      
      tab.cont.i.j.k.ref <- subset(tab.cont.i.j.ref, cat.AI == k)
      
      if(dim(tab.cont.i.j.k.ref)[1] == 0){
        names1 <- names(tab.cont.i.j.k.ref)
        line1 <- c("1970_2000","historical", i, j, k, 0)
        tab.cont.i.j.k.ref <- rbind(tab.cont.i.j.k.ref, line1) %>%
          setNames(names1)
        
      }else{
        tab.cont.i.j.k.ref <- tab.cont.i.j.k.ref
      }
      
      tab.cont.i.j.k.ref <- tab.cont.i.j.k.ref %>% rbind(tab.cont.i.j.k.ref[rep(1,lines),])
      tab.cont.i.j.k$diff <- tab.cont.i.j.k$count - as.numeric(tab.cont.i.j.k.ref$count)  
      
      list_ref[[index]] <- tab.cont.i.j.k
    }
  }
}


tab.cont.diff <- bind_rows(list_ref, .id = "column_label") %>%
  group_by(period, model, Continent, cat.AI) %>% 
  summarise(diff.mean = mean(diff, na.rm = T), diff.sd = sd(diff, na.rm = T))
write.table(tab.cont.diff, "tab.cont.diff.txt")

tab.cont.diff <- read.table("tab.cont.diff.txt")

tab_list_percent <- list()
tab_list_sd <- list()

for(i in unique(tab.cont.diff$Continent)){
  tab.i <- subset(tab.cont.diff, Continent == i)
  
  df.percent <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "diff.mean")
  
  df.sd <- tab.i %>% reshape2::dcast(period + model ~ cat.AI, value.var = "diff.sd")
  # sd for a sum: square root of sum of squared sd
  
  drycats <- which(names(df.percent) %in% c("Hyperarid","Arid","Semi-arid","Dry subhumid"))
  
  df.percent <- df.percent %>% mutate("Sum drylands" = rowSums(.[drycats], na.rm = T))
  df.sd <- df.sd %>% mutate("Sum drylands" = sqrt(rowSums((.[drycats])^2, na.rm = T)),2)
  
  missing <- setdiff(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"), names(df.percent))
  
  df.percent[, missing] <- 0
  df.sd[,missing] <- 0
  
  df.percent <- df.percent %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.percent[is.na(df.percent)] <- 0
  df.percent[,c(3:9)] <- round(df.percent[,c(3:9)], 1)
  
  tab_list_percent[[i]] <- df.percent
  
  df.sd <- df.sd %>% select(c("model","period","Hyperarid","Arid","Semi-arid","Dry subhumid","Sum drylands","Humid","Cold","NA"))
  df.sd[is.na(df.sd)] <- 0
  df.sd[ , c(3:9)] <- round(df.sd[, c(3:9)], 1)
  
  tab_list_sd[[i]] <- df.sd
}

k_list <- list()

for(i in names(tab_list_percent)){
  
  df <- tab_list_percent[[i]] %>% select(c("model", "period"))
  
  df$Hyperarid <- paste(tab_list_percent[[i]]$Hyperarid, tab_list_sd[[i]]$Hyperarid, sep = " ± ")
  df$Arid <- paste(tab_list_percent[[i]]$Arid, tab_list_sd[[i]]$Arid, sep = " ± ")
  df$'Semi-Arid' <- paste(tab_list_percent[[i]][,5], tab_list_sd[[i]][,5], sep = " ± ")
  df$'Dry subhumid' <- paste(tab_list_percent[[i]][,6], tab_list_sd[[i]][,6], sep = " ± ")
  df$'Sum drylands' <- paste(tab_list_percent[[i]][,7], tab_list_sd[[i]][,7], sep = " ± ")
  df$Humid <- paste(tab_list_percent[[i]]$Humid, tab_list_sd[[i]]$Humid, sep = " ± ")
  df$Cold <- paste(tab_list_percent[[i]]$Cold, tab_list_sd[[i]]$Cold, sep = " ± ")
  
  k_list[[i]] <- kable(df[order(df$model, decreasing = F)], caption = i, digits = 2) %>% kable_styling(bootstrap_options = "bordered") %>%
    column_spec(8, italic = T, background = colorvec[5]) %>% row_spec(2, bold = T) 
  
}

k_list

$AFRICA

AFRICA

model	period	Hyperarid	Arid	Sum drylands	Semi-Arid	Humid	Dry subhumid	Cold
historical	1850_1880	-5.8 ± 24.1	3 ± 22.7	1.5 ± 48.9	4.8 ± 34	2.8 ± 38.7	-0.4 ± 12	0 ± 0
historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
historical	1985_2015	-10.5 ± 19.5	11 ± 21.6	2.8 ± 33.1	5.9 ± 11.9	-2.7 ± 15.5	-3.5 ± 10.1	-2 ± 0
SSP245	2030_2060	-10.5 ± 41	24.3 ± 29.7	42.2 ± 66.6	25.8 ± 40.3	-37.8 ± 53.2	2.5 ± 15.5	0 ± 0
SSP370	2030_2060	-29.5 ± 88.4	-5.2 ± 60.9	23.7 ± 145	45.5 ± 86.1	-21.4 ± 43.4	12.9 ± 45.6	-30 ± 0
SSP585	2030_2060	-7.8 ± 50.9	13.6 ± 39.6	39 ± 82.2	33.2 ± 47.6	-34.6 ± 52.7	0 ± 18.4	0 ± 0
SSP245	2070_2100	10.1 ± 66.8	9 ± 46.2	58.8 ± 99	42.1 ± 51.9	-54.5 ± 61.7	-2.3 ± 22.7	0 ± 0
SSP370	2070_2100	19.5 ± 66.3	18.5 ± 58.2	69.2 ± 103.5	31.1 ± 46.6	-67.1 ± 52.3	0.2 ± 27.4	-28 ± 0
SSP585	2070_2100	45.6 ± 209.6	9.2 ± 79.8	67.5 ± 240.3	16.5 ± 83.8	-62.8 ± 98	-3.8 ± 20.8	0 ± 0

$ASIA

ASIA

model	period	Hyperarid	Arid	Sum drylands	Semi-Arid	Humid	Dry subhumid	Cold
historical	1850_1880	4.3 ± 28	8.8 ± 41.5	40.2 ± 88.6	11.3 ± 61.5	90.2 ± 176.8	15.8 ± 39.6	-130.5 ± 287.7
historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
historical	1985_2015	-5.2 ± 20	-12.6 ± 42	-36.4 ± 57.8	-17 ± 29.8	-33.6 ± 169.1	-1.6 ± 17.2	70 ± 234.1
SSP245	2030_2060	-15.7 ± 38.1	24.2 ± 40.5	57.8 ± 123.7	22.2 ± 88.9	192.7 ± 202.9	27 ± 65.5	-250.5 ± 375.8
SSP370	2030_2060	-23.8 ± 61.4	21.6 ± 45.7	-4.5 ± 98.3	0 ± 55.9	174.2 ± 120.1	-2.3 ± 26	-169.7 ± 100.4
SSP585	2030_2060	-1.8 ± 30.8	23.8 ± 34.1	72.9 ± 115.9	27.1 ± 91.2	250 ± 201.3	23.8 ± 54.9	-322.9 ± 358.7
SSP245	2070_2100	-8.5 ± 31.9	23.6 ± 37.9	68.1 ± 133.8	30.1 ± 101.2	294.5 ± 208.7	22.8 ± 72.2	-362.5 ± 381.8
SSP370	2070_2100	40.1 ± 64.3	25.3 ± 63.8	101.2 ± 108.9	32.7 ± 53	241.2 ± 179.9	3.1 ± 29	-342.3 ± 156
SSP585	2070_2100	55.1 ± 178.4	59.5 ± 146	215 ± 313.2	77.9 ± 202.8	324.1 ± 262.6	22.5 ± 61.9	-539.1 ± 416.4

$CENTRAL-AMERICA

CENTRAL-AMERICA

model	period	Hyperarid	Arid	Sum drylands	Semi-Arid	Humid	Dry subhumid	Cold
historical	1850_1880	0 ± 0	-0.3 ± 4.7	2.5 ± 21.5	-0.5 ± 18	-2 ± 28.2	3.3 ± 10.8	0 ± 0
historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
historical	1985_2015	0 ± 0	1.7 ± 4	4.2 ± 8.6	2.5 ± 6.3	-4 ± 5	-0.1 ± 4.3	-2 ± 0
SSP245	2030_2060	0 ± 0	9.8 ± 11.1	24 ± 29.7	17.3 ± 25.7	-23.5 ± 26.8	-3.2 ± 10	0 ± 0
SSP370	2030_2060	0 ± 0	11.8 ± 12.4	15.1 ± 31.2	3.3 ± 24	-14.6 ± 37.3	0 ± 15.6	0 ± 0
SSP585	2030_2060	0 ± 0	12.8 ± 14	31.5 ± 36.2	21.1 ± 29.1	-31 ± 25.5	-2.5 ± 16.4	0 ± 0
SSP245	2070_2100	0 ± 0	11 ± 10.3	31.9 ± 36.4	22.6 ± 29.2	-31.5 ± 27.3	-1.7 ± 19	0 ± 0
SSP370	2070_2100	12 ± 15.6	23.2 ± 15.1	50.2 ± 39	18.7 ± 29.5	-39.6 ± 23.8	-3.7 ± 13.4	0 ± 0
SSP585	2070_2100	6 ± 7.1	25.2 ± 24.9	60.6 ± 42.8	35.6 ± 28.6	-55.1 ± 42.4	-6.2 ± 18.6	0 ± 0

$EUROPE

EUROPE

model	period	Hyperarid	Arid	Sum drylands	Semi-Arid	Humid	Dry subhumid	Cold
historical	1850_1880	0 ± 0	4.7 ± 11.2	20.7 ± 20.7	8.2 ± 13.3	50.3 ± 106.6	7.8 ± 11.3	-70 ± 115
historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
historical	1985_2015	0 ± 0	-2.2 ± 9.8	-9.2 ± 27.8	-5.8 ± 24.6	12.9 ± 48.5	-1.2 ± 8.5	-4.3 ± 73.1
SSP245	2030_2060	0 ± 0	4.7 ± 7.4	33.7 ± 18.9	17.5 ± 14.6	97.5 ± 126	11.5 ± 9.6	-128.5 ± 146.3
SSP370	2030_2060	0 ± 0	7.2 ± 8.5	35.2 ± 16.1	17.6 ± 11.2	90.5 ± 75.1	10.4 ± 7.8	-125.2 ± 73.5
SSP585	2030_2060	0 ± 0	9.8 ± 8.2	41.1 ± 18.4	19.1 ± 14.2	125.9 ± 134.4	12.2 ± 8.4	-164.8 ± 144.4
SSP245	2070_2100	0 ± 0	8.6 ± 11.2	44.2 ± 24.7	21.6 ± 19.9	123.1 ± 120.8	14 ± 9.3	-166.6 ± 139.5
SSP370	2070_2100	0 ± 0	19.2 ± 13.7	86.9 ± 44.6	37.7 ± 28	92.2 ± 110	30 ± 31.8	-177.7 ± 82.2
SSP585	2070_2100	0 ± 0	18.3 ± 19	81.8 ± 39.7	37.4 ± 24.7	150.7 ± 131	26.1 ± 24.6	-230.9 ± 148.4

$EUROPE-AFRICA

EUROPE-AFRICA

model	period	Hyperarid	Arid	Sum drylands	Semi-Arid	Humid	Dry subhumid	Cold
historical	1850_1880	-2.1 ± 5.2	-1.2 ± 9	-1.8 ± 19.9	1.8 ± 14.9	1 ± 20.7	-0.4 ± 8.2	13 ± 0
historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
historical	1985_2015	0.5 ± 5.8	0.2 ± 8.2	-1.4 ± 24.4	-4.2 ± 17.9	2.5 ± 16.7	2 ± 13.2	-6.5 ± 9.2
SSP245	2030_2060	11.2 ± 8.1	7 ± 10	26.5 ± 25.8	11.6 ± 19.8	-22.2 ± 35.5	-3.3 ± 10.1	-54 ± 0
SSP370	2030_2060	14 ± 6.8	6.1 ± 9.9	26.1 ± 20.3	10.5 ± 14.3	-23.1 ± 20.7	-4.5 ± 8	-56 ± 0
SSP585	2030_2060	16.9 ± 9.9	7.8 ± 10.2	35.5 ± 25.8	16.5 ± 17.2	-33.2 ± 32	-5.8 ± 13	-57 ± 0
SSP245	2070_2100	19.4 ± 8.2	6.8 ± 11.4	38.6 ± 29	14.9 ± 23.1	-36.5 ± 35.9	-2.5 ± 10.7	0 ± 0
SSP370	2070_2100	25.3 ± 16.7	-0.9 ± 23.1	35.5 ± 41.9	17.7 ± 22.1	-30.9 ± 28.5	-6.5 ± 21.3	0 ± 0
SSP585	2070_2100	32.4 ± 12.3	11.3 ± 19.9	56.3 ± 36.9	21.2 ± 22.1	-55.5 ± 27.4	-8.7 ± 18.1	0 ± 0

$NORTH-AMERICA

NORTH-AMERICA

model	period	Hyperarid	Arid	Sum drylands	Semi-Arid	Humid	Dry subhumid	Cold
historical	1850_1880	0 ± 0	0.4 ± 3.1	37 ± 71.9	31.8 ± 66.3	21.7 ± 91.3	4.9 ± 27.7	-56.1 ± 174
historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
historical	1985_2015	0 ± 0	0.3 ± 4.2	-5.4 ± 48.8	-3 ± 32.2	-19.8 ± 70	-2.7 ± 36.4	25.1 ± 114
SSP245	2030_2060	0 ± 0	3.4 ± 2.9	77.5 ± 109.8	60 ± 104.2	82.2 ± 94.5	14.2 ± 34.7	-159.2 ± 214.3
SSP370	2030_2060	0 ± 0	5.5 ± 6.9	62 ± 56.5	44.5 ± 45.2	70.5 ± 48.5	12.1 ± 33.2	-131.6 ± 77
SSP585	2030_2060	0 ± 0	7.1 ± 7.9	99 ± 102.3	62.1 ± 91	103.4 ± 97.9	29.8 ± 46.2	-201.3 ± 203.2
SSP245	2070_2100	0 ± 0	4.7 ± 3.1	87.1 ± 99.3	56.9 ± 86	145.9 ± 122.4	25.5 ± 49.7	-231.9 ± 215.6
SSP370	2070_2100	4 ± 0	15.7 ± 20	139.6 ± 80.2	91.5 ± 59.6	108 ± 110.3	28.4 ± 49.8	-243.9 ± 106
SSP585	2070_2100	5 ± 0	20.2 ± 38.4	147.1 ± 145.7	91.2 ± 136	198.8 ± 99.8	30.8 ± 35.4	-339.6 ± 220.2

$OCEANIA

OCEANIA

model	period	Hyperarid	Arid	Sum drylands	Semi-Arid	Humid	Dry subhumid	Cold
historical	1850_1880	1 ± 0	6 ± 42.8	4.5 ± 57.7	-1.2 ± 38.1	-3 ± 12.1	-1.2 ± 6.8	0.2 ± 1.2
historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
historical	1985_2015	1 ± 0	5.1 ± 32.9	5.6 ± 40.3	-2.2 ± 21.3	-3 ± 7.9	1.8 ± 9.4	-0.2 ± 1
SSP245	2030_2060	45 ± 0	50.2 ± 62.1	62.8 ± 74.9	-26.5 ± 39.9	-19.8 ± 18.3	-5.8 ± 12.4	-1.5 ± 0.7
SSP370	2030_2060	47 ± 0	42.4 ± 50.3	71.7 ± 61.2	-12.4 ± 33.3	-27 ± 61.2	-5.3 ± 10.3	-3.3 ± 3.2
SSP585	2030_2060	23.5 ± 0.7	33.8 ± 55.9	33.3 ± 71.6	-23.2 ± 43.1	-12 ± 20.1	-0.8 ± 12.2	-1.5 ± 0.7
SSP245	2070_2100	45 ± 0	66.6 ± 81.5	60.7 ± 101.2	-48.5 ± 58.5	-17.4 ± 18.8	-2.5 ± 13	-5 ± 0
SSP370	2070_2100	70.7 ± 84	63.2 ± 110.8	70 ± 153.2	-57 ± 62.9	-14.2 ± 57.5	-6.8 ± 13.5	-5 ± 0
SSP585	2070_2100	60.3 ± 48.7	69.8 ± 126.8	71.2 ± 164.5	-50 ± 91.8	-23.9 ± 22.1	-9 ± 14.2	0 ± 0

$SOUTH-AMERICA

SOUTH-AMERICA

model	period	Hyperarid	Arid	Sum drylands	Semi-Arid	Humid	Dry subhumid	Cold
historical	1850_1880	-0.1 ± 0.4	5.1 ± 13.9	26.3 ± 86.6	20.1 ± 83.4	-20.6 ± 79.5	1.3 ± 18.8	-6.8 ± 35.1
historical	1970_2000	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0	0 ± 0
historical	1985_2015	0.5 ± 1.2	2.7 ± 5.9	16.7 ± 19.2	10.5 ± 15.7	-14.2 ± 13.2	3 ± 9.3	-2.5 ± 5.8
SSP245	2030_2060	1.6 ± 1.1	13.2 ± 26.6	93.8 ± 100	55.6 ± 93.1	-78.5 ± 97.2	23.5 ± 24.8	-18.8 ± 42
SSP370	2030_2060	1.2 ± 0.9	5.9 ± 19.7	69.4 ± 52.4	38.2 ± 31.9	-59.5 ± 55.6	24.1 ± 36.6	-12.9 ± 18.7
SSP585	2030_2060	1.5 ± 1.1	19.2 ± 29.9	123.7 ± 97.4	72.8 ± 84.5	-107.5 ± 99.4	30.2 ± 38.1	-20.1 ± 41.2
SSP245	2070_2100	1.4 ± 1	22.6 ± 49.9	124 ± 111.5	68.9 ± 90.8	-107 ± 130.1	31 ± 41.3	-22.8 ± 43.7
SSP370	2070_2100	2 ± 2.3	22.9 ± 23.8	121.2 ± 92.7	73.5 ± 82.1	-109.5 ± 109.3	22.7 ± 35.9	-15 ± 15.1
SSP585	2070_2100	7.4 ± 15.3	40.5 ± 51.6	182.8 ± 140.1	96.9 ± 122.8	-160.2 ± 133.1	38 ± 40.7	-29.4 ± 47.8

10 Changes in proportion of aridity categories, figure

10.1 Mean changes of sum drylands by continent, period and model

list_plots <- list()

for(i in names(tab_list_percent)){
  df.percent <- tab_list_percent[[i]]
  df.sd <- tab_list_sd[[i]]
  
  df <- merge(df.percent, df.sd, by = c("model", "period"))
  
  g <- ggplot(data = df)+geom_point(aes(x=period,y= get("Sum drylands.x"), col = model), shape = "\u2605", size = 5, position = position_dodge(width = 0.5))+
    geom_errorbar(aes(x = period, ymin = get("Sum drylands.x")-get("Sum drylands.y"), ymax = get("Sum drylands.x")+get("Sum drylands.y"), col = model), position = position_dodge(width = 0.5))+
    scale_color_manual(values = colorvec[c(1,3,9,11)])+
    ylim(0,100)+
    labs(x="", y = "Percent of drylands", title = i)+
    theme_minimal()
  list_plots[[i]] <- g
}

ggarrange(plotlist = list_plots, ncol = 2, nrow = 4, common.legend = T, legend = "bottom")

10.2 Colplot by continent

tab.flow <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 ) %>%
  group_by(period, model, Continent, cat.AI) %>%
  summarise(count = n()) %>%
  ungroup() %>% group_by(period, model, Continent) %>% mutate(percent = round(count/sum(count)*100, 1), count = NULL)

df245 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP245")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b245 <- ggplot(data = df245, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  facet_grid(rows = vars(Continent), switch = "y")+
  scale_y_continuous(position = "right")+
  labs(title = "SSP 2-4.5", y = "%", x = "")+
  theme_minimal()

df370 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP370")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b370 <- ggplot(data = df370, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  facet_grid(rows = vars(Continent), switch = "y")+
  scale_y_continuous(position = "right")+  labs(title = "SSP 3-7.0", y = "%", x = "")+
  theme_minimal()


df585 <- rbind(subset(tab.flow, model == "historical" & period %in% c("1850_1880", "1970_2000")), subset(tab.flow, model == "SSP585")) %>% mutate(lab.y = (rev(cumsum(rev(percent)))) - percent*0.5) %>% 
  subset(cat.AI %in% c("Hyperarid", "Arid", "Semi-arid", "Dry subhumid", "Humid","Cold"))

b585 <- ggplot(data = df585, aes(x = period, y = percent))+
  geom_col(aes(group = cat.AI, col = cat.AI, fill = cat.AI))+
  geom_label(aes(y = lab.y, label = percent))+
  scale_color_manual(values = col.cat, aesthetics = c("col", "fill"), na.translate = F)+
  facet_grid(rows = vars(Continent), switch = "y")+
  scale_y_continuous(position = "right")+  labs(title = "SSP 5-8.5", y = "%", x = "")+
  theme_minimal()


ggarrange(plotlist = list(b245, b370, b585), common.legend = T, legend = "bottom", ncol = 3)

10.3 Evolution of aridity index by Continent and period

tab.percent.region <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period != "1985_2015") %>%
  group_by(period, model, Continent, Name) %>% select(c("diff.AI.mean.percent")) %>%
  summarise(AI.diff.mean = mean(diff.AI.mean.percent, na.rm = T), AI.diff.sd= sd(diff.AI.mean.percent, na.rm = T)) %>%
  ungroup()

ipcc_regions.sf <- st_as_sf(ipcc_regions)

percent.region.sf <- merge(ipcc_regions.sf, tab.percent.region, all.y = T)

pbreaks <- c(-10, 0, 10, 20, 30, 40)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "azure","cornsilk", "#f4a582")
#colscale <- c("#c6dbef","white","#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#993344")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot(data = subset(percent.region.sf, period == j & model == i)) + geom_sf(aes(fill = AI.diff.mean))+
          binned_scale(aesthetics = "fill", breaks = pbreaks, palette = function(x) rev(colscale))+
      borders(colour = "grey60")+
      geom_sf_label(aes(label = round(AI.diff.mean,1)),size = 2) +
      labs(title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

for(i in unique(percent.region.sf$Continent)){

  df <- percent.region.sf %>% filter(Continent == i & period != "1850_1880") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

  dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  
    
g <- ggplot(dfref)+geom_point(aes(x = period, y = AI.diff.mean, col = Name))+
  geom_line(aes(x = period, y = AI.diff.mean, col = Name, group = Name))+
  scale_color_manual(values = c(colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef","#e8ecff", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")))+
  facet_grid(cols = vars(model), scales = "free_x")+
  theme_minimal()+
  theme(legend.position = "right")
  
plot(g)
}

tab.percent.continent <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period != "1985_2015") %>%
  group_by(period, model, Continent) %>% select(c("diff.AI.mean.percent")) %>%
  summarise(AI.diff.mean = mean(diff.AI.mean.percent, na.rm = T), AI.diff.sd= sd(diff.AI.mean.percent, na.rm = T)) %>%
  ungroup()


df <- tab.percent.continent %>% filter(period != "1850_1880" & Continent != "POLAR") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  

dfref$Continent <- str_replace(dfref$Continent, "-", "\n")
    
ggplot(dfref, aes(x = period, y = AI.diff.mean))+
  geom_point(aes(size = AI.diff.sd, col = Continent))+
  geom_line(aes(group = Continent), size = 0.1)+
  facet_grid(rows = vars(model), cols = vars(Continent))+
  scale_size_continuous(breaks = c(0,1,5,10,100), transform = "pseudo_log", name = "Standard deviation")+
  scale_color_viridis_d(end = 0.8)+
  labs(y = "Difference in aridity index compared to 1970-2000, %", x = "")+
  guides(color = F)+
  ylim(-40,25)+
  theme_minimal()+
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 90))

tab.percent.continent <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period != "1985_2015") %>%
  group_by(period, model, Name) %>% select(c("diff.AI.mean.percent")) %>%
  summarise(AI.diff.mean = mean(diff.AI.mean.percent, na.rm = T), AI.diff.sd= sd(diff.AI.mean.percent, na.rm = T)) %>%
  ungroup()


df <- tab.percent.continent %>% filter(period != "1850_1880") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  

dfref$Name <- str_replace(dfref$Name, "-", "\n") 

ggplot(dfref, aes(x = period, y = AI.diff.mean))+
  geom_point(aes(size = AI.diff.sd, col = Name))+
  geom_line(aes(group = Name), size = 0.1)+
  facet_grid(cols = vars(model), rows = vars(Name), switch = "both")+
  scale_size_continuous(breaks = c(0,1,5,10,100), transform = "pseudo_log", name = "Standard deviation")+
  scale_color_viridis_d()+
labs(y = "Difference in aridity index compared to 1970-2000, %", x = "")+
  guides(color = F)+
  theme_minimal()+
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 90))+  
  scale_y_continuous(position = "right", limits = c(-40,30))

tab.percent.continent <- cmip6s %>% 
  subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period %in% c("1970_2000","2030_2060","2070_2100"))

df <- tab.percent.continent %>% filter(Continent != "POLAR") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  

dfref$Continent <- str_replace(dfref$Continent, "-", "\n")

ggplot(dfref, aes(x = period, y = diff.AI.mean.percent))+
  geom_boxplot(aes(col = Continent), outliers = F)+
  facet_grid(rows = vars(model), cols = vars(Continent), space = "free_y")+
  scale_color_viridis_d(end = 0.8)+
  labs(y = "Difference in aridity index compared to 1970-2000, %", x = "")+
  guides(color = F)+
  theme_minimal()+
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 90))

cmip6s$pr.percent <- with(cmip6s, diff.pr.mean/spr.ref.mean*100)

tab.percent.continent <- cmip6s %>% subset(!Continent %in% c("SOUTHERN","PACIFIC","ATLANTIC","INDIAN","ARCTIC") & lat > -55 & period != "1985_2015") %>%
  group_by(period, model, Continent) %>% select(c("diff.t.mean", "pr.percent")) %>%
  summarise(t.diff.mean = mean(diff.t.mean, na.rm = T), t.diff.sd= sd(diff.t.mean, na.rm = T),
            pr.diff.mean = mean(pr.percent, na.rm = T), pr.diff.sd= sd(pr.percent, na.rm = T)) %>%
  ungroup()


df <- tab.percent.continent %>% filter(period != "1850_1880" & Continent != "POLAR") 
ref <- df[which(df$model == "historical"),]
ref1 <- ref %>% mutate(model = "SSP245")
ref2 <- ref %>% mutate(model = "SSP370")
ref3 <- ref %>% mutate(model = "SSP585")

dfref <- rbind(subset(df, model != "historical"), ref1) %>% rbind(ref2) %>% rbind(ref3)  

dfref$Continent <- str_replace(dfref$Continent, "-", "\n")
    
ggplot(dfref, aes(x = period))+
  geom_point(aes(y = t.diff.mean*3, size = t.diff.sd, col = "Temperature, °C"), alpha = 0.5)+
  geom_line(aes(y = t.diff.mean*3, group = Continent, col = "Temperature, °C"), size = 0.1)+
  
  geom_point(aes(y = pr.diff.mean, size = pr.diff.sd, col = "Precipitation, %"), alpha = 0.5)+
  geom_line(aes(y = pr.diff.mean, group = Continent, col = "Precipitation, %"), size = 0.1)+
  
  scale_color_manual(values = c("#3182bd", "#d6604d"))+
  facet_grid(rows = vars(model), cols = vars(Continent))+
  labs(y = "Precipitation anomalies, %", x = "", size = "Standard deviation", color = "")+
  ylim(-45,25)+
  theme_minimal()+
  theme(legend.position = "bottom", axis.text.x = element_text(angle = 90))+
  
  scale_y_continuous(sec.axis = sec_axis(trans = ~./3, name = "Temperature anomalies, °C"))

11 Changes of variables over time and scenarios

11.1 Temperature

11.1.1 Absolute temperature

quantile(cmip6s$t.mean, probs = seq(0,1,0.1)) 
##          0%         10%         20%         30%         40%         50% 
## -53.4587313 -42.5828261 -28.4644864 -17.7502317  -7.2080554  -0.2580962 
##         60%         70%         80%         90%        100% 
##   7.1198406  15.3318920  22.8072968  26.4971129  34.3727997

t.breaks <- c(-50,-40,-30,-20,-10,0,10,20,30)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")

map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = t.mean))+
    borders(colour = "grey60")+
    binned_scale(aesthetics = "fill", breaks = t.breaks, palette = function(x) colscale,
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "Mean annual temperature, °C")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = t.mean))+
      borders(colour = "grey60")+
      binned_scale(aesthetics = "fill", breaks = t.breaks, palette = function(x) colscale,
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Mean annual temperature, °C", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}



ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

#labels = c("2030-2060", "2070-2100"

11.1.2 Temperature anomalies

quantile(cmip6s$diff.t.mean, probs = seq(0,1,0.1)) 
##         0%        10%        20%        30%        40%        50%        60% 
## -3.5775507 -0.5616741  0.0000000  1.1670755  1.8544034  2.2566573  2.6600679 
##        70%        80%        90%       100% 
##  3.1444949  3.7454555  4.5875606 10.7627615

t.breaks <- c(-3,-0.5,0,1.5,2,3,4,5,10)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")

map_list <- list()

for(i in c("1850_1880", "1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = diff.t.mean))+
    borders(colour = "grey60")+
   scale_fill_gradient2(low = "#08519c",mid = "white", high = "#b2182b", midpoint = 0, transform = "pseudo_log", breaks = c(-2, -1, 0, 1, 2, 3, 4))+
      labs(fill = "Precipitation anomalies in %", title = i)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_tile(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = diff.t.mean))+
      borders(colour = "grey60")+
   scale_fill_gradient2(low = "#08519c",mid = "white", high = "#b2182b", midpoint = 0, transform = "pseudo_log", breaks = c(-2, -1, 0, 1, 2, 3, 4))+
      labs(fill = "Temperature anomalies in °C", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

### Binned anomalies

YOR <- brewer.pal(n = 11, name ="YlOrRd")
coltemp <- c("[0,1]" = YOR[1], "(1,2]" = YOR[2], "(2,3]" = YOR[3], "(3,4]" = YOR[4], "(4,5]" = YOR[5], "(5,6]" = YOR[6], "(6,7]" = YOR[7], "(7,8]" = YOR[8], "(8,9]" = YOR[9], "(9,10]" = "grey20")


map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot(data = subset(cmip6s, period == j & model == i), 
                aes(x=lon, y = lat,  fill = cut(diff.t.mean, breaks = c(0,1,2,3,4,5,6,7,8,9,10), include.lowest = T))) + 
      geom_tile(show.legend = T)+
      borders(colour = "grey60")+
      scale_fill_manual(values = coltemp, drop = F)+
      labs(fill = "Temperature\nanomalies in °C", title = index)+
      theme_void()+ylim(-55,90)+
     theme(legend.position = "bottom", legend.text =  element_text(angle = 0))+
    guides(fill = guide_legend(nrow = 2))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

11.2 Precipitation

11.2.1 Absolute values of precipitations

quantile(cmip6s$pr.mean, probs = seq(0,1,0.1)) 
##          0%         10%         20%         30%         40%         50% 
##    8.579696   55.000319  104.132310  200.852004  337.042710  463.946055 
##         60%         70%         80%         90%        100% 
##  579.741473  712.099280  959.490392 1431.315602 6863.723388

p.breaks <- c(10,50,100,200,300,400,500,600,700)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")

map_list <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = pr.mean))+
    borders(colour = "grey60")+
    binned_scale(aesthetics = "fill", breaks = p.breaks, palette = function(x) rev(colscale),
                 guide = guide_legend(label.theme = element_text(angle = 0)))+
    labs(title = i, fill = "Annual mean\nprecipitation, mm")+
    theme_void()+ylim(-55,90)+
    theme(legend.position = "right")
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = pr.mean))+
      borders(colour = "grey60")+
      binned_scale(aesthetics = "fill", breaks = p.breaks, palette = function(x) rev(colscale),
                   guide = guide_legend(label.theme = element_text(angle = 0)))+
      labs(fill = "Annual mean\nprecipitation, mm", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom")
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

11.2.2 Precipitation anomalies

cmip6s$pr.percent <- with(cmip6s, diff.pr.mean/spr.ref.mean*100)

quantile(cmip6s$pr.percent, probs = seq(0,1,0.1)) 
##           0%          10%          20%          30%          40%          50% 
## -42.02999493  -2.87352237   0.00000000   0.04663897   2.46136429   5.19906939 
##          60%          70%          80%          90%         100% 
##   9.34748067  13.67239987  18.11758968  24.02085190 148.15226168

p.breaks <- c(-43,-3,0,3,5,9,13,18,24,148)
colscale <- c("#08519c", "#3182bd", "#6baed6", "#9ecae1", "#c6dbef", "#fddbc7", "#f4a582", "#d6604d", "#b2182b", "#67001f")

map_list <- list()

for(i in c("1850_1880","1970_2000","1985_2015")){
  g <- ggplot() + geom_raster(data = subset(cmip6s, period == i & model == "historical"), aes(x=lon, y = lat,  fill = pr.percent))+
    borders(colour = "grey60")+
   scale_fill_gradient2(high = "#08519c", mid = "white", low = "#b2182b", midpoint = 0, transform = "pseudo_log", breaks = c(-10, -5, -1, 0, 1, 5, 10, 15, 30))+
      labs(fill = "Precipitation anomalies in %", title = i)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
  map_list[[i]] <- g
}

ggpubr::ggarrange(plotlist = map_list, ncol = 1, nrow = 3, common.legend = T, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = pr.percent))+
      borders(colour = "grey60")+
   scale_fill_gradient2(high = "#08519c", mid = "white", low = "#b2182b", midpoint = 0, transform = "pseudo_log", breaks = c(-10, -5, -1, 0, 1, 5, 10, 15, 30))+
      labs(fill = "Precipitation anomalies in %", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")

map_list <- list()

for(i in c("SSP245","SSP370","SSP585")){
  for (j in c("2030_2060","2070_2100")){
    index <- paste(i, j, sep = " , ")
    g <- ggplot() + geom_raster(data = subset(cmip6s, period == j & model == i), aes(x=lon, y = lat,  fill = pr.percent))+
      borders(colour = "grey60")+
   scale_fill_binned_divergingx(mid = 0, palette = "RdBu", rev = F, n_interp = 7, breaks = c(-20, -10, 0, 10, 20, 30, 40, 50, 60))+
      labs(fill = "Precipitation anomalies in %", title = index)+
      theme_void()+ylim(-55,90)+
      theme(legend.position = "bottom", legend.text =  element_text(angle = 0))
    map_list[[index]] <- g
  }}

ggpubr::ggarrange(plotlist = map_list,ncol = 2, nrow = 3, common.legend = T, legend = "bottom")